diff --git a/dumux/assembly/staggeredlocalassembler.hh b/dumux/assembly/staggeredlocalassembler.hh
index 8351df9199bbc05c29ea33aefcd64170591fa160..a725718f766d574b8cbaaee0fd2ce6dfb43f27ea 100644
--- a/dumux/assembly/staggeredlocalassembler.hh
+++ b/dumux/assembly/staggeredlocalassembler.hh
@@ -78,6 +78,9 @@ class StaggeredLocalAssembler<TypeTag,
     using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
     using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
     using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
+    using FaceVariables = typename GET_PROP_TYPE(TypeTag, FaceVariables);
+    using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
 
     static constexpr bool enableGlobalFluxVarsCache = GET_PROP_VALUE(TypeTag, EnableGlobalFluxVariablesCache);
 
@@ -111,8 +114,10 @@ public:
         auto elemFluxVarsCache = localView(gridVariables.gridFluxVarsCache());
         elemFluxVarsCache.bind(element, fvGeometry, curElemVolVars);
 
-        auto& curGlobalFaceVars = gridVariables.curGridFaceVars();
-        auto& prevGlobalFaceVars = gridVariables.prevGridFaceVars();
+        auto curElemeFaceVars = localView(gridVariables.curGridFaceVars());
+        auto prevElemeFaceVars = localView(gridVariables.prevGridFaceVars());
+        curElemeFaceVars.bind(element, fvGeometry, curSol);
+        prevElemeFaceVars.bind(element, fvGeometry, curSol); //TODO: stationary
 
         const bool isStationary = localResidual.isStationary();
         auto prevElemVolVars = localView(gridVariables.prevGridVolVars());
@@ -128,8 +133,8 @@ public:
                                                                              fvGeometry,
                                                                              prevElemVolVars,
                                                                              curElemVolVars,
-                                                                             prevGlobalFaceVars,
-                                                                             curGlobalFaceVars,
+                                                                             prevElemeFaceVars,
+                                                                             curElemeFaceVars,
                                                                              elemBcTypes,
                                                                              elemFluxVarsCache)[0];
 
@@ -148,8 +153,8 @@ public:
                                                                             scvf,
                                                                             prevElemVolVars,
                                                                             curElemVolVars,
-                                                                            prevGlobalFaceVars,
-                                                                            curGlobalFaceVars,
+                                                                            prevElemeFaceVars,
+                                                                            curElemeFaceVars,
                                                                             elemBcTypes,
                                                                             elemFluxVarsCache)[0];
 
@@ -163,8 +168,8 @@ public:
                                 curSol,
                                 prevElemVolVars,
                                 curElemVolVars,
-                                prevGlobalFaceVars,
-                                curGlobalFaceVars,
+                                prevElemeFaceVars,
+                                curElemeFaceVars,
                                 elemFluxVarsCache,
                                 elemBcTypes,
                                 jac,
@@ -261,8 +266,8 @@ protected:
                                         const SolutionVector& curSol,
                                         const ElementVolumeVariables& prevElemVolVars,
                                         ElementVolumeVariables& curElemVolVars,
-                                        const GlobalFaceVars& prevGlobalFaceVars,
-                                        GlobalFaceVars& curGlobalFaceVars,
+                                        const ElementFaceVariables& prevElemFaceVars,
+                                        ElementFaceVariables& curElemFaceVars,
                                         ElementFluxVariablesCache& elemFluxVarsCache,
                                         const ElementBoundaryTypes& elemBcTypes,
                                         JacobianMatrix& matrix,
@@ -270,16 +275,16 @@ protected:
                                         const FaceSolutionVector& faceResidualCache)
 {
     // compute the derivatives of the cell center residuals with respect to cell center dofs
-    dCCdCC_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, ccResidual);
+    dCCdCC_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, elemFluxVarsCache, elemBcTypes, matrix, ccResidual);
 
     // compute the derivatives of the cell center residuals with respect to face dofs
-    dCCdFace_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, ccResidual);
+    dCCdFace_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, elemFluxVarsCache, elemBcTypes, matrix, ccResidual);
 
     // compute the derivatives of the face residuals with respect to cell center dofs
-    dFacedCC_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, faceResidualCache);
+    dFacedCC_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, elemFluxVarsCache, elemBcTypes, matrix, faceResidualCache);
 
     // compute the derivatives of the face residuals with respect to face dofs
-    dFacedFace_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevGlobalFaceVars, curGlobalFaceVars, elemFluxVarsCache, elemBcTypes, matrix, faceResidualCache);
+    dFacedFace_(assembler, element, fvGeometry, curSol, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, elemFluxVarsCache, elemBcTypes, matrix, faceResidualCache);
 }
 
 /*!
@@ -292,8 +297,8 @@ static void dCCdCC_(Assembler& assembler,
              const SolutionVector& curSol,
              const ElementVolumeVariables& prevElemVolVars,
              ElementVolumeVariables& curElemVolVars,
-             const GlobalFaceVars& prevGlobalFaceVars,
-             const GlobalFaceVars& curGlobalFaceVars,
+             const ElementFaceVariables& prevElemFaceVars,
+             const ElementFaceVariables& curElemFaceVars,
              ElementFluxVariablesCache& elemFluxVarsCache,
              const ElementBoundaryTypes& elemBcTypes,
              JacobianMatrix& matrix,
@@ -330,7 +335,7 @@ static void dCCdCC_(Assembler& assembler,
            curVolVars.update(elemSol, problem, elementJ, scvJ);
 
           const auto deflectedResidual = localResidual.evalCellCenter(problem, element, fvGeometry, prevElemVolVars, curElemVolVars,
-                                         prevGlobalFaceVars, curGlobalFaceVars,
+                                         prevElemFaceVars, curElemFaceVars,
                                          elemBcTypes, elemFluxVarsCache);
 
            const auto partialDeriv = (deflectedResidual - ccResidual) / eps;
@@ -354,8 +359,8 @@ static void dCCdFace_(Assembler& assembler,
                       const SolutionVector& curSol,
                       const ElementVolumeVariables& prevElemVolVars,
                       const ElementVolumeVariables& curElemVolVars,
-                      const GlobalFaceVars& prevGlobalFaceVars,
-                      GlobalFaceVars& curGlobalFaceVars,
+                      const ElementFaceVariables& prevElemFaceVars,
+                      ElementFaceVariables& curElemFaceVars,
                       ElementFluxVariablesCache& elemFluxVarsCache,
                       const ElementBoundaryTypes& elemBcTypes,
                       JacobianMatrix& matrix,
@@ -368,6 +373,7 @@ static void dCCdFace_(Assembler& assembler,
 
    const auto& problem = assembler.problem();
    auto& localResidual = assembler.localResidual();
+   auto& gridVariables = assembler.gridVariables();
 
    // build derivatives with for cell center dofs w.r.t. cell center dofs
    const auto cellCenterGlobalI = assembler.fvGridGeometry().elementMapper().index(element);
@@ -377,8 +383,8 @@ static void dCCdFace_(Assembler& assembler,
         const auto globalJ = scvfJ.dofIndex();
 
        // get the faceVars of the face with respect to which we are going to build the derivative
-       auto origFaceVars = curGlobalFaceVars.faceVars(scvfJ.index());
-       auto& curFaceVars = curGlobalFaceVars.faceVars(scvfJ.index());
+       auto& faceVars = getFaceVarAccess(gridVariables.curGridFaceVars(), curElemFaceVars, scvfJ);
+       const auto origFaceVars(faceVars);
 
        for(auto pvIdx : PriVarIndices(faceIdx))
        {
@@ -391,11 +397,11 @@ static void dCCdFace_(Assembler& assembler,
            priVars[pvIdx] += eps;
 
            faceSolution[globalJ][pvIdx] += eps;
-           curFaceVars.update(scvfJ, faceSolution);
+           faceVars.update(scvfJ, faceSolution);
 
            const auto deflectedResidual = localResidual.evalCellCenter(problem, element, fvGeometry,
                                           prevElemVolVars, curElemVolVars,
-                                          prevGlobalFaceVars, curGlobalFaceVars,
+                                          prevElemFaceVars, curElemFaceVars,
                                           elemBcTypes, elemFluxVarsCache);
 
            const auto partialDeriv = (deflectedResidual - ccResidual) / eps;
@@ -404,7 +410,7 @@ static void dCCdFace_(Assembler& assembler,
            updateGlobalJacobian_(matrix[cellCenterIdx][faceIdx], cellCenterGlobalI, globalJ, pvIdx - Indices::faceOffset, partialDeriv);
 
            // restore the original faceVars
-           curFaceVars = origFaceVars;
+           faceVars = origFaceVars;
        }
    }
 }
@@ -419,8 +425,8 @@ static void dFacedCC_(Assembler& assembler,
                       const SolutionVector& curSol,
                       const ElementVolumeVariables& prevElemVolVars,
                       ElementVolumeVariables& curElemVolVars,
-                      const GlobalFaceVars& prevGlobalFaceVars,
-                      const GlobalFaceVars& curGlobalFaceVars,
+                      const ElementFaceVariables& prevElemFaceVars,
+                      const ElementFaceVariables& curElemFaceVars,
                       ElementFluxVariablesCache& elemFluxVarsCache,
                       const ElementBoundaryTypes& elemBcTypes,
                       JacobianMatrix& matrix,
@@ -461,7 +467,7 @@ static void dFacedCC_(Assembler& assembler,
 
                const auto deflectedResidual = localResidual.evalFace(problem, element, fvGeometry, scvf,
                                               prevElemVolVars, curElemVolVars,
-                                              prevGlobalFaceVars, curGlobalFaceVars,
+                                              prevElemFaceVars, curElemFaceVars,
                                               elemBcTypes, elemFluxVarsCache);
 
                const auto partialDeriv = (deflectedResidual - cachedResidual[scvf.localFaceIdx()]) / eps;
@@ -486,8 +492,8 @@ static void dFacedFace_(Assembler& assembler,
                         const SolutionVector& curSol,
                         const ElementVolumeVariables& prevElemVolVars,
                         const ElementVolumeVariables& curElemVolVars,
-                        const GlobalFaceVars& prevGlobalFaceVars,
-                        GlobalFaceVars& curGlobalFaceVars,
+                        const ElementFaceVariables& prevElemFaceVars,
+                        ElementFaceVariables& curElemFaceVars,
                         ElementFluxVariablesCache& elemFluxVarsCache,
                         const ElementBoundaryTypes& elemBcTypes,
                         JacobianMatrix& matrix,
@@ -499,6 +505,7 @@ static void dFacedFace_(Assembler& assembler,
     const auto& problem = assembler.problem();
     auto& localResidual = assembler.localResidual();
     const auto& connectivityMap = assembler.fvGridGeometry().connectivityMap();
+    auto& gridVariables = assembler.gridVariables();
 
    for(auto&& scvf : scvfs(fvGeometry))
    {
@@ -509,8 +516,8 @@ static void dFacedFace_(Assembler& assembler,
        for(const auto& globalJ : connectivityMap(faceIdx, faceIdx, scvf.index()))
        {
            // get the faceVars of the face with respect to which we are going to build the derivative
-           auto origFaceVars = curGlobalFaceVars.faceVars(scvf.index());
-           auto& curFaceVars = curGlobalFaceVars.faceVars(scvf.index());
+        auto& faceVars = getFaceVarAccess(gridVariables.curGridFaceVars(), curElemFaceVars, scvf);
+        const auto origFaceVars(faceVars);
 
            for(auto pvIdx : PriVarIndices(faceIdx))
            {
@@ -519,11 +526,11 @@ static void dFacedFace_(Assembler& assembler,
                const Scalar eps = numericEpsilon(faceSolution[globalJ][pvIdx], faceIdx, faceIdx);
 
                faceSolution[globalJ][pvIdx] += eps;
-               curFaceVars.update(scvf, faceSolution);
+               faceVars.update(scvf, faceSolution);
 
                const auto deflectedResidual = localResidual.evalFace(problem, element, fvGeometry, scvf,
                                               prevElemVolVars, curElemVolVars,
-                                              prevGlobalFaceVars, curGlobalFaceVars,
+                                              prevElemFaceVars, curElemFaceVars,
                                               elemBcTypes, elemFluxVarsCache);
 
                const auto partialDeriv = (deflectedResidual - cachedResidual[scvf.localFaceIdx()]) / eps;
@@ -532,7 +539,7 @@ static void dFacedFace_(Assembler& assembler,
                updateGlobalJacobian_(matrix[faceIdx][faceIdx], faceGlobalI, globalJ, pvIdx - Indices::faceOffset, partialDeriv);
 
                // restore the original faceVars
-               curFaceVars = origFaceVars;
+               faceVars = origFaceVars;
            }
        }
    }
@@ -603,6 +610,16 @@ private:
     static typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), VolumeVariables&>::type
     getVolVarAccess(GridVolumeVariables& gridVolVars, ElementVolumeVariables& elemVolVars, const SubControlVolume& scv)
     { return gridVolVars.volVars(scv); }
+
+    template<class T = TypeTag>
+    static typename std::enable_if<!GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), FaceVariables&>::type //TODO: introduce correct property
+    getFaceVarAccess(GlobalFaceVars& gridFaceVars, ElementFaceVariables& elemFaceVars, const SubControlVolumeFace& scvf)
+    { return elemFaceVars[scvf]; }
+
+    template<class T = TypeTag>
+    static typename std::enable_if<GET_PROP_VALUE(T, EnableGlobalVolumeVariablesCache), FaceVariables&>::type
+    getFaceVarAccess(GlobalFaceVars& gridFaceVars, ElementFaceVariables& elemFaceVars, const SubControlVolumeFace& scvf)
+    { return gridFaceVars.faceVars(scvf.index()); }
 };
 
 } // end namespace Dumux
diff --git a/dumux/discretization/staggered/elementfacevariables.hh b/dumux/discretization/staggered/elementfacevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ee142f93e237e562734f46de50f27927d23c2844
--- /dev/null
+++ b/dumux/discretization/staggered/elementfacevariables.hh
@@ -0,0 +1,76 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ * \brief The face variables class for free flow staggered grid models
+ */
+#ifndef DUMUX_DISCRETIZATION_STAGGERED_ELEMENTFACEVARIABLES_HH
+#define DUMUX_DISCRETIZATION_STAGGERED_ELEMENTFACEVARIABLES_HH
+
+#include <dumux/common/basicproperties.hh>
+
+namespace Dumux
+{
+
+
+template<class TypeTag>
+class StaggeredElementFaceVariables
+{
+    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+    using Element = typename GridView::template Codim<0>::Entity;
+    using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
+    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+    using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
+    using FaceVariables = typename GET_PROP_TYPE(TypeTag, FaceVariables);
+    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
+    using IndexType = typename GridView::IndexSet::IndexType;
+
+public:
+
+    StaggeredElementFaceVariables(const GlobalFaceVars& globalFacesVars) : globalFaceVarsPtr_(&globalFacesVars) {}
+
+    const FaceVariables& operator [](const SubControlVolumeFace& scvf) const
+    { return globalFaceVars().faceVars(scvf.index()); }
+
+    // operator for the access with an index
+    // needed for cc methods for the access to the boundary volume variables
+    const FaceVariables& operator [](const IndexType scvfIdx) const
+    { return globalFaceVars().faceVars(scvfIdx); }
+
+
+    //! For compatibility reasons with the case of not storing the vol vars.
+    //! function to be called before assembling an element, preparing the vol vars within the stencil
+    void bind(const Element& element,
+              const FVElementGeometry& fvGeometry,
+              const SolutionVector& sol)
+    {}
+
+
+    //! The global volume variables object we are a restriction of
+    const GlobalFaceVars& globalFaceVars() const
+    { return *globalFaceVarsPtr_; }
+
+
+private:
+    const GlobalFaceVars* globalFaceVarsPtr_;
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/discretization/staggered/freeflow/facevariables.hh b/dumux/discretization/staggered/freeflow/facevariables.hh
index e8bdf8b5d34e4755ea8c22b81f0bb198e06d9709..a4e80fcc88ef95c95fb93930993e652ae4476739 100644
--- a/dumux/discretization/staggered/freeflow/facevariables.hh
+++ b/dumux/discretization/staggered/freeflow/facevariables.hh
@@ -40,6 +40,7 @@ class StaggeredFaceVariables
     using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
     using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
     using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
+    using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
     using FaceSolution = typename GET_PROP_TYPE(TypeTag, StaggeredFaceSolution);
 
     struct SubFaceData
@@ -51,6 +52,19 @@ class StaggeredFaceVariables
     };
 
 public:
+    StaggeredFaceVariables() = default;
+
+    StaggeredFaceVariables(const GlobalFaceVars& globalFacesVars) : globalFaceVarsPtr_(&globalFacesVars) {}
+
+    const StaggeredFaceVariables& operator [](const SubControlVolumeFace& scvf) const
+    { return globalFaceVars().faceVars(scvf.index()); }
+
+    // // operator for the access with an index
+    // // needed for cc methods for the access to the boundary volume variables
+    // const StaggeredFaceVariables& operator [](const IndexType scvIdx) const
+    // { return globalVolVars().volVars(scvIdx); }
+
+
     void update(const FacePrimaryVariables &facePrivars)
     {
         velocitySelf_ = facePrivars;
@@ -120,8 +134,13 @@ public:
         return subFaceVelocities_[localSubFaceIdx];
     }
 
+    //! The global volume variables object we are a restriction of
+    const GlobalFaceVars& globalFaceVars() const
+    { return *globalFaceVarsPtr_; }
+
 
 private:
+    const GlobalFaceVars* globalFaceVarsPtr_;
     Scalar velocity_;
     Scalar velocitySelf_;
     Scalar velocityOpposite_;
diff --git a/dumux/discretization/staggered/globalfacevariables.hh b/dumux/discretization/staggered/globalfacevariables.hh
index 0c0caa72342648530bdab0ce708c2742f67393a1..3f8edb66db81e367e72621e29fa37ca495303446 100644
--- a/dumux/discretization/staggered/globalfacevariables.hh
+++ b/dumux/discretization/staggered/globalfacevariables.hh
@@ -29,6 +29,11 @@
 namespace Dumux
 {
 
+namespace Properties
+{
+    NEW_PROP_TAG(ElementFaceVariables);
+}
+
 template<class TypeTag>
 class StaggeredGlobalFaceVariables
 {
@@ -37,6 +42,7 @@ class StaggeredGlobalFaceVariables
     using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
     using FaceVariables = typename GET_PROP_TYPE(TypeTag, FaceVariables);
+    using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
     using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry);
     using IndexType = typename GridView::IndexSet::IndexType;
 
@@ -49,14 +55,7 @@ public:
 
     void update(const FVGridGeometry& fvGridGeometry, const SolutionVector& sol)
     {
-
         const auto& faceSol = sol[faceIdx];
-
-
-
-        // const auto numFaceDofs = fvGridGeometry.gridView().size(1);
-
-        // faceVariables_.resize(numFaceDofs);
         faceVariables_.resize(fvGridGeometry.numScvf());
 
         for(auto&& element : elements(fvGridGeometry.gridView()))
@@ -67,21 +66,8 @@ public:
             for(auto&& scvf : scvfs(fvGeometry))
             {
                 faceVariables_[scvf.index()].update(scvf, faceSol);
-                auto test = StaggeredFaceSolution<TypeTag>(scvf, faceSol, fvGridGeometry);
-                // std::cout << "test" << test[scvf.dofIndex()] << std::endl;
             }
         }
-
-        // for(auto& faceVariable : faceVariables_)
-        // {
-        //     faceVariable.update()
-        // }
-        // // assert(faceVariables_.size() == faceSol.size());
-        //
-        // for(int i = 0; i < numFaceDofs; ++i)
-        // {
-        //     faceVariables_[i].update(faceSol[i]);
-        // }
     }
 
     const FaceVariables& faceVars(const IndexType facetIdx) const
@@ -89,6 +75,16 @@ public:
 
     FaceVariables& faceVars(const IndexType facetIdx)
     { return faceVariables_[facetIdx]; }
+
+
+    /*!
+     * \brief Return a local restriction of this global object
+     *        The local object is only functional after calling its bind/bindElement method
+     *        This is a free function that will be found by means of ADL
+     */
+    friend inline ElementFaceVariables localView(const StaggeredGlobalFaceVariables& global)
+    { return ElementFaceVariables(global); }
+
 private:
     const Problem& problem_() const
     { return *problemPtr_; }
diff --git a/dumux/discretization/staggered/properties.hh b/dumux/discretization/staggered/properties.hh
index 67cb5b475cbb0e703565de36bf828824605d7ea4..15f153ebeeaca4399846da48fa7866ab3ff9979f 100644
--- a/dumux/discretization/staggered/properties.hh
+++ b/dumux/discretization/staggered/properties.hh
@@ -45,6 +45,7 @@
 #include <dumux/discretization/staggered/fvelementgeometry.hh>
 #include <dumux/discretization/staggered/globalfacevariables.hh>
 #include <dumux/discretization/staggered/facesolution.hh>
+#include <dumux/discretization/staggered/elementfacevariables.hh>
 
 #include <dumux/common/intersectionmapper.hh>
 #include <dune/istl/multitypeblockvector.hh>
@@ -62,6 +63,7 @@ namespace Properties
 NEW_PROP_TAG(CellCenterSolutionVector);
 NEW_PROP_TAG(FaceSolutionVector);
 NEW_PROP_TAG(StaggeredFaceSolution);
+NEW_PROP_TAG(ElementFaceVariables);
 
 //! Type tag for the box scheme.
 NEW_TYPE_TAG(StaggeredModel, INHERITS_FROM(FiniteVolumeModel, NumericModel));
@@ -113,6 +115,8 @@ SET_TYPE_PROP(StaggeredModel, IntersectionMapper, Dumux::ConformingGridIntersect
 
 SET_TYPE_PROP(StaggeredModel, StaggeredFaceSolution, StaggeredFaceSolution<TypeTag>);
 
+SET_TYPE_PROP(StaggeredModel, ElementFaceVariables, StaggeredElementFaceVariables<TypeTag>);
+
 //! Definition of the indices for cell center and face dofs in the global solution vector
 SET_PROP(StaggeredModel, DofTypeIndices)
 {
diff --git a/dumux/freeflow/staggered/fluxvariables.hh b/dumux/freeflow/staggered/fluxvariables.hh
index a3de0a7ccc9a3f44fb9ba2a336963468653ce758..a19d1220c26b92dfe45ddf839ef39bb250d504a4 100644
--- a/dumux/freeflow/staggered/fluxvariables.hh
+++ b/dumux/freeflow/staggered/fluxvariables.hh
@@ -66,7 +66,6 @@ class FreeFlowFluxVariablesImpl<TypeTag, false>
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
-    using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
     using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
     using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
     using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
@@ -74,6 +73,8 @@ class FreeFlowFluxVariablesImpl<TypeTag, false>
     using IndexType = typename GridView::IndexSet::IndexType;
     using Stencil = std::vector<IndexType>;
 
+    using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
+
     static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms);
 
     enum {
@@ -100,7 +101,7 @@ public:
                                                         const Element &element,
                                                         const FVElementGeometry& fvGeometry,
                                                         const ElementVolumeVariables& elemVolVars,
-                                                        const GlobalFaceVars& globalFaceVars,
+                                                        const ElementFaceVariables& elemFaceVars,
                                                         const SubControlVolumeFace &scvf,
                                                         const FluxVariablesCache& fluxVarsCache)
     {
@@ -111,8 +112,7 @@ public:
         const auto& outsideVolVars = scvf.boundary() ?  insideVolVars : elemVolVars[scvf.outsideScvIdx()];
 
         CellCenterPrimaryVariables flux(0.0);
-        // const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
-        const Scalar velocity = globalFaceVars.faceVars(scvf.index()).velocitySelf();
+        const Scalar velocity = elemFaceVars[scvf].velocitySelf();
 
         const bool insideIsUpstream = sign(scvf.outerNormalScalar()) == sign(velocity) ? true : false;
         const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
@@ -191,19 +191,19 @@ public:
     * \param scvf The sub control volume face
     * \param fvGeometry The finite-volume geometry
     * \param elemVolVars All volume variables for the element
-    * \param globalFaceVars The face variables
+    * \param elementFaceVars The face variables
     */
    FacePrimaryVariables computeNormalMomentumFlux(const Problem& problem,
                                                   const Element& element,
                                                   const SubControlVolumeFace& scvf,
                                                   const FVElementGeometry& fvGeometry,
                                                   const ElementVolumeVariables& elemVolVars,
-                                                  const GlobalFaceVars& globalFaceVars)
+                                                  const ElementFaceVariables& elementFaceVars)
    {
        const auto insideScvIdx = scvf.insideScvIdx();
        const auto& insideVolVars = elemVolVars[insideScvIdx];
-       const Scalar velocitySelf = globalFaceVars.faceVars(scvf.index()).velocitySelf() ;
-       const Scalar velocityOpposite = globalFaceVars.faceVars(scvf.index()).velocityOpposite();
+       const Scalar velocitySelf = elementFaceVars[scvf].velocitySelf() ;
+       const Scalar velocityOpposite = elementFaceVars[scvf].velocityOpposite();
        FacePrimaryVariables normalFlux(0.0);
 
        if(navierStokes)
@@ -243,28 +243,20 @@ public:
    * \param scvf The sub control volume face
    * \param fvGeometry The finite-volume geometry
    * \param elemVolVars All volume variables for the element
-   * \param globalFaceVars The face variables
+   * \param elementFaceVars The face variables
    */
   FacePrimaryVariables computeTangetialMomentumFlux(const Problem& problem,
                                                     const Element& element,
                                                     const SubControlVolumeFace& scvf,
                                                     const FVElementGeometry& fvGeometry,
                                                     const ElementVolumeVariables& elemVolVars,
-                                                    const GlobalFaceVars& globalFaceVars)
+                                                    const ElementFaceVariables& elementFaceVars)
   {
       FacePrimaryVariables tangentialFlux(0.0);
-
-    //   // convenience function to get the velocity on a face
-    //   auto velocity = [&globalFaceVars](const int dofIdx)
-    //   {
-    //       return globalFaceVars.faceVars(scvf.dofIdx()).velocity();
-    //   };
-    auto& faceVars = globalFaceVars.faceVars(scvf.index());
-
-    const int numSubFaces = scvf.pairData().size();
+      auto& faceVars = elementFaceVars[scvf];
+      const int numSubFaces = scvf.pairData().size();
 
       // account for all sub-faces
-    //   for(const auto& subFaceData : scvf.pairData())
       for(int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx)
       {
           const auto eIdx = scvf.insideScvIdx();
@@ -305,7 +297,6 @@ private:
                                                                      const FaceVars& faceVars,
                                                                      const int localSubFaceIdx)
   {
-    //   const Scalar transportingVelocity = velocity(subFaceData.normalPair.first);
       const Scalar transportingVelocity = faceVars.subFaceData(localSubFaceIdx).velocityNormalInside;
       const auto insideScvIdx = normalFace.insideScvIdx();
       const auto outsideScvIdx = normalFace.outsideScvIdx();
@@ -324,12 +315,10 @@ private:
 
       if(innerElementIsUpstream)
           transportedVelocity = faceVars.subFaceData(localSubFaceIdx).velocityParallelInside;
-        //   transportedVelocity = velocity(scvf.dofIndex());
       else
       {
           const int outerDofIdx = scvf.pairData(localSubFaceIdx).outerParallelFaceDofIdx;
           if(outerDofIdx >= 0)
-            //   transportedVelocity = velocity(outerDofIdx);
               transportedVelocity = faceVars.subFaceData(localSubFaceIdx).velocityParallelOutside;
           else // this is the case when the outer parallal dof would lie outside the domain TODO: discuss which one is better
             //   transportedVelocity = problem.dirichlet(makeGhostFace(subFaceData.virtualOuterParallelFaceDofPos))[faceIdx][scvf.directionIndex()];
@@ -370,10 +359,8 @@ private:
       const Scalar muAvg = (insideVolVars.viscosity() + outsideVolVars.viscosity()) * 0.5;
 
       // the normal derivative
-    //   const int innerNormalVelocityIdx = subFaceData.normalPair.first;
       const int outerNormalVelocityIdx = scvf.pairData(localSubFaceIdx).normalPair.second;
 
-    //   const Scalar innerNormalVelocity = velocity(innerNormalVelocityIdx);
       const Scalar innerNormalVelocity = faceVars.subFaceData(localSubFaceIdx).velocityNormalInside;
 
       const Scalar outerNormalVelocity = outerNormalVelocityIdx >= 0 ?
diff --git a/dumux/freeflow/staggered/localresidual.hh b/dumux/freeflow/staggered/localresidual.hh
index 606f5fe2874d91a195840c8b5018c89761ac56f2..e4557974b06b65a57058c662dfa656ebed584c12 100644
--- a/dumux/freeflow/staggered/localresidual.hh
+++ b/dumux/freeflow/staggered/localresidual.hh
@@ -35,6 +35,7 @@ NEW_PROP_TAG(EnableInertiaTerms);
 NEW_PROP_TAG(ReplaceCompEqIdx);
 NEW_PROP_TAG(EnergyFluxVariables);
 NEW_PROP_TAG(NormalizePressure);
+NEW_PROP_TAG(ElementFaceVariables);
 }
 
 /*!
@@ -86,6 +87,7 @@ class StaggeredNavierStokesResidualImpl<TypeTag, false> : public Dumux::Staggere
     using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
     using EnergyLocalResidual = typename GET_PROP_TYPE(TypeTag, EnergyLocalResidual);
     using EnergyFluxVariables = typename GET_PROP_TYPE(TypeTag, EnergyFluxVariables);
+    using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
 
 
     using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
@@ -109,7 +111,6 @@ class StaggeredNavierStokesResidualImpl<TypeTag, false> : public Dumux::Staggere
     };
 
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
-    using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
 
     static constexpr bool navierStokes = GET_PROP_VALUE(TypeTag, EnableInertiaTerms);
     static constexpr bool normalizePressure = GET_PROP_VALUE(TypeTag, NormalizePressure);
@@ -123,15 +124,15 @@ public:
                                                         const Element &element,
                                                         const FVElementGeometry& fvGeometry,
                                                         const ElementVolumeVariables& elemVolVars,
-                                                        const GlobalFaceVars& globalFaceVars,
+                                                        const ElementFaceVariables& elemFaceVars,
                                                         const SubControlVolumeFace &scvf,
                                                         const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
         FluxVariables fluxVars;
         CellCenterPrimaryVariables flux = fluxVars.computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars,
-                                                 globalFaceVars, scvf, elemFluxVarsCache[scvf]);
+                                                 elemFaceVars, scvf, elemFluxVarsCache[scvf]);
 
-        EnergyFluxVariables::energyFlux(flux, problem, element, fvGeometry, elemVolVars, globalFaceVars, scvf, elemFluxVarsCache[scvf]);
+        EnergyFluxVariables::energyFlux(flux, problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache[scvf]);
 
         return flux;
     }
@@ -140,7 +141,7 @@ public:
                                                           const Element &element,
                                                           const FVElementGeometry& fvGeometry,
                                                           const ElementVolumeVariables& elemVolVars,
-                                                          const GlobalFaceVars& globalFaceVars,
+                                                          const ElementFaceVariables& elemFaceVars,
                                                           const SubControlVolume &scv) const
     {
         return CellCenterPrimaryVariables(0.0);
@@ -179,10 +180,10 @@ public:
      */
     FacePrimaryVariables computeStorageForFace(const SubControlVolumeFace& scvf,
                                                const VolumeVariables& volVars,
-                                               const GlobalFaceVars& globalFaceVars)
+                                               const ElementFaceVariables& elementFaceVars)
     {
         FacePrimaryVariables storage(0.0);
-        const Scalar velocity = globalFaceVars.faceVars(scvf.index()).velocity();
+        const Scalar velocity = elementFaceVars[scvf].velocity();
         storage[0] = volVars.density() * velocity;
         return storage;
     }
@@ -190,7 +191,7 @@ public:
     FacePrimaryVariables computeSourceForFace(const Problem& problem,
                                               const SubControlVolumeFace& scvf,
                                               const ElementVolumeVariables& elemVolVars,
-                                              const GlobalFaceVars& globalFaceVars) const
+                                              const ElementFaceVariables& elementFaceVars) const
     {
         FacePrimaryVariables source(0.0);
         const auto insideScvIdx = scvf.insideScvIdx();
@@ -207,21 +208,21 @@ public:
      * \param scvf The sub control volume face
      * \param fvGeometry The finite-volume geometry
      * \param elemVolVars All volume variables for the element
-     * \param globalFaceVars The face variables
+     * \param elementFaceVars The face variables
      */
     FacePrimaryVariables computeFluxForFace(const Problem& problem,
                                             const Element& element,
                                             const SubControlVolumeFace& scvf,
                                             const FVElementGeometry& fvGeometry,
                                             const ElementVolumeVariables& elemVolVars,
-                                            const GlobalFaceVars& globalFaceVars,
+                                            const ElementFaceVariables& elementFaceVars,
                                             const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
         FacePrimaryVariables flux(0.0);
         FluxVariables fluxVars;
-        flux += fluxVars.computeNormalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, globalFaceVars);
-        flux += fluxVars.computeTangetialMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, globalFaceVars);
-        flux += computePressureTerm_(problem, element, scvf, fvGeometry, elemVolVars, globalFaceVars);
+        flux += fluxVars.computeNormalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elementFaceVars);
+        flux += fluxVars.computeTangetialMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elementFaceVars);
+        flux += computePressureTerm_(problem, element, scvf, fvGeometry, elemVolVars, elementFaceVars);
         return flux;
     }
 
@@ -233,14 +234,14 @@ protected:
     void evalBoundary_(const Element& element,
                        const FVElementGeometry& fvGeometry,
                        const ElementVolumeVariables& elemVolVars,
-                       const GlobalFaceVars& faceVars,
+                       const ElementFaceVariables& elemFaceVars,
                        const ElementBoundaryTypes& elemBcTypes,
                        const ElementFluxVariablesCache& elemFluxVarsCache)
     {
-        evalBoundaryForCellCenter_(element, fvGeometry, elemVolVars, faceVars, elemBcTypes, elemFluxVarsCache);
+        evalBoundaryForCellCenter_(element, fvGeometry, elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache);
         for (auto&& scvf : scvfs(fvGeometry))
         {
-            evalBoundaryForFace_(element, fvGeometry, scvf, elemVolVars, faceVars, elemBcTypes, elemFluxVarsCache);
+            evalBoundaryForFace_(element, fvGeometry, scvf, elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache);
         }
     }
 
@@ -252,7 +253,7 @@ protected:
                                     const Element& element,
                                     const FVElementGeometry& fvGeometry,
                                     const ElementVolumeVariables& elemVolVars,
-                                    const GlobalFaceVars& faceVars,
+                                    const ElementFaceVariables& elemFaceVars,
                                     const ElementBoundaryTypes& elemBcTypes,
                                     const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
@@ -260,7 +261,7 @@ protected:
         {
             if (scvf.boundary())
             {
-                auto boundaryFlux = computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, faceVars, scvf, elemFluxVarsCache);
+                auto boundaryFlux = computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);
 
                 // handle the actual boundary conditions:
                 const auto bcTypes = problem.boundaryTypes(element, scvf);
@@ -315,7 +316,7 @@ protected:
                               const FVElementGeometry& fvGeometry,
                               const SubControlVolumeFace& scvf,
                               const ElementVolumeVariables& elemVolVars,
-                              const GlobalFaceVars& faceVars,
+                              const ElementFaceVariables& elementFaceVars,
                               const ElementBoundaryTypes& elemBcTypes,
                               const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
@@ -328,7 +329,7 @@ protected:
             if(bcTypes.isDirichlet(momentumBalanceIdx))
             {
                 // const Scalar velocity = faceVars.faceVars(scvf.dofIndex()).velocity();
-                const Scalar velocity = faceVars.faceVars(scvf.index()).velocitySelf();
+                const Scalar velocity = elementFaceVars[scvf].velocitySelf();
                 const Scalar dirichletValue = problem.dirichlet(element, scvf)[faceIdx][scvf.directionIndex()];
                 residual = velocity - dirichletValue;
             }
@@ -338,7 +339,7 @@ protected:
             if(bcTypes.isSymmetry())
             {
                 // const Scalar velocity = faceVars.faceVars(scvf.dofIndex()).velocity();
-                const Scalar velocity = faceVars.faceVars(scvf.index()).velocitySelf();
+                const Scalar velocity = elementFaceVars[scvf].velocitySelf();
                 const Scalar fixedValue = 0.0;
                 residual = velocity - fixedValue;
             }
@@ -347,7 +348,7 @@ protected:
             if(bcTypes.isOutflow(momentumBalanceIdx))
             {
                 if(bcTypes.isDirichlet(massBalanceIdx))
-                    residual += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, faceVars, elemFluxVarsCache);
+                    residual += computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elementFaceVars, elemFluxVarsCache);
                 else
                     DUNE_THROW(Dune::InvalidStateException, "Face at " << scvf.center()  << " has an outflow BC for the momentum balance but no Dirichlet BC for the pressure!");
             }
@@ -364,14 +365,14 @@ private:
      * \param scvf The sub control volume face
      * \param fvGeometry The finite-volume geometry
      * \param elemVolVars All volume variables for the element
-     * \param globalFaceVars The face variables
+     * \param elementFaceVars The face variables
      */
     FacePrimaryVariables computePressureTerm_(const Problem& problem,
                                               const Element& element,
                                               const SubControlVolumeFace& scvf,
                                               const FVElementGeometry& fvGeometry,
                                               const ElementVolumeVariables& elemVolVars,
-                                              const GlobalFaceVars& globalFaceVars) const
+                                              const ElementFaceVariables& elementFaceVars) const
     {
         const auto insideScvIdx = scvf.insideScvIdx();
         const auto& insideVolVars = elemVolVars[insideScvIdx];
diff --git a/dumux/freeflow/staggerednc/fluxvariables.hh b/dumux/freeflow/staggerednc/fluxvariables.hh
index ae4a1c6574fc62824519dc7137bf5c5fd55f838c..83542c6eb03eda7bd065b43733ea38d6a7cac072 100644
--- a/dumux/freeflow/staggerednc/fluxvariables.hh
+++ b/dumux/freeflow/staggerednc/fluxvariables.hh
@@ -36,6 +36,7 @@ namespace Properties
 NEW_PROP_TAG(EnableComponentTransport);
 NEW_PROP_TAG(EnableEnergyBalance);
 NEW_PROP_TAG(EnableInertiaTerms);
+NEW_PROP_TAG(ElementFaceVariables);
 }
 
 // // forward declaration
@@ -60,7 +61,7 @@ class FreeFlowFluxVariablesImpl<TypeTag, true> : public FreeFlowFluxVariablesImp
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
-    using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
+    using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
     using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
     using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
     using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
@@ -99,13 +100,13 @@ public:
                                                         const Element &element,
                                                         const FVElementGeometry& fvGeometry,
                                                         const ElementVolumeVariables& elemVolVars,
-                                                        const GlobalFaceVars& globalFaceVars,
+                                                        const ElementFaceVariables& elemFaceVars,
                                                         const SubControlVolumeFace &scvf,
                                                         const FluxVariablesCache& fluxVarsCache)
     {
         CellCenterPrimaryVariables flux(0.0);
 
-        flux += advectiveFluxForCellCenter_(problem, fvGeometry, elemVolVars, globalFaceVars, scvf);
+        flux += advectiveFluxForCellCenter_(problem, fvGeometry, elemVolVars, elemFaceVars, scvf);
         flux += MolecularDiffusionType::diffusiveFluxForCellCenter(problem, fvGeometry, elemVolVars, scvf);
         return flux;
     }
@@ -115,7 +116,7 @@ private:
     CellCenterPrimaryVariables advectiveFluxForCellCenter_(const Problem& problem,
                                                           const FVElementGeometry& fvGeometry,
                                                           const ElementVolumeVariables& elemVolVars,
-                                                          const GlobalFaceVars& globalFaceVars,
+                                                          const ElementFaceVariables& elemFaceVars,
                                                           const SubControlVolumeFace &scvf)
     {
         CellCenterPrimaryVariables flux(0.0);
@@ -123,7 +124,7 @@ private:
         const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
         const auto& insideVolVars = elemVolVars[insideScv];
 
-        const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
+        const Scalar velocity = elemFaceVars[scvf].velocitySelf();
 
         const bool insideIsUpstream = sign(scvf.outerNormalScalar()) == sign(velocity);
         const Scalar upWindWeight = GET_PROP_VALUE(TypeTag, ImplicitUpwindWeight);
diff --git a/dumux/freeflow/staggerednc/localresidual.hh b/dumux/freeflow/staggerednc/localresidual.hh
index 0fb4538ede8936bb1c70cdedda1fc4d08b8e96ca..db16c348a81be626c2daa7da8f7d996d625f571a 100644
--- a/dumux/freeflow/staggerednc/localresidual.hh
+++ b/dumux/freeflow/staggerednc/localresidual.hh
@@ -87,7 +87,6 @@ class StaggeredNavierStokesResidualImpl<TypeTag, true> : public StaggeredNavierS
     };
 
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
-    using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
     using EnergyLocalResidual = typename GET_PROP_TYPE(TypeTag, EnergyLocalResidual);
 
     static constexpr int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
diff --git a/dumux/freeflow/staggeredni/fluxvariables.hh b/dumux/freeflow/staggeredni/fluxvariables.hh
index 8fa221d1a522ff7e8f8e1931832d4c52c3c5785f..028f37129637114fd6186e50f69c2c60c4f32350 100644
--- a/dumux/freeflow/staggeredni/fluxvariables.hh
+++ b/dumux/freeflow/staggeredni/fluxvariables.hh
@@ -28,6 +28,11 @@
 namespace Dumux
 {
 
+namespace Properties
+{
+    NEW_PROP_TAG(ElementFaceVariables);
+}
+
 /*!
  * \ingroup ImplicitModel
  * \brief The flux variables class
@@ -52,7 +57,7 @@ class FreeFlowEnergyFluxVariablesImplementation<TypeTag, false>
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
-    using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
+    using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
     using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
     using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
     using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
@@ -64,7 +69,7 @@ public:
                            const Element &element,
                            const FVElementGeometry& fvGeometry,
                            const ElementVolumeVariables& elemVolVars,
-                           const GlobalFaceVars& globalFaceVars,
+                           const ElementFaceVariables& elemFaceVars,
                            const SubControlVolumeFace &scvf,
                            const FluxVariablesCache& fluxVarsCache)
     { }
@@ -82,7 +87,7 @@ class FreeFlowEnergyFluxVariablesImplementation<TypeTag, true>
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
-    using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
+    using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
     using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
     using FluxVariablesCache = typename GET_PROP_TYPE(TypeTag, FluxVariablesCache);
     using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
@@ -98,11 +103,11 @@ public:
                            const Element &element,
                            const FVElementGeometry& fvGeometry,
                            const ElementVolumeVariables& elemVolVars,
-                           const GlobalFaceVars& globalFaceVars,
+                           const ElementFaceVariables& elemFaceVars,
                            const SubControlVolumeFace &scvf,
                            const FluxVariablesCache& fluxVarsCache)
     {
-        flux[energyBalanceIdx] += advectiveFluxForCellCenter_(problem, fvGeometry, elemVolVars, globalFaceVars, scvf);
+        flux[energyBalanceIdx] += advectiveFluxForCellCenter_(problem, fvGeometry, elemVolVars, elemFaceVars, scvf);
         flux[energyBalanceIdx] += HeatConductionType::diffusiveFluxForCellCenter(problem, element, fvGeometry, elemVolVars, scvf);
     }
 
@@ -112,13 +117,13 @@ private:
     * \param problem The problem
     * \param fvGeometry The finite-volume geometry
     * \param elemVolVars All volume variables for the element
-    * \param globalFaceVars The face variables
+    * \param elemFaceVars The face variables
     * \param scvf The sub control volume face
     */
     static Scalar advectiveFluxForCellCenter_(const Problem& problem,
                                               const FVElementGeometry& fvGeometry,
                                               const ElementVolumeVariables& elemVolVars,
-                                              const GlobalFaceVars& globalFaceVars,
+                                              const ElementFaceVariables& elemFaceVars,
                                               const SubControlVolumeFace &scvf)
     {
         Scalar flux(0.0);
@@ -138,7 +143,7 @@ private:
 
         const auto& outsideVolVars = isOutflow ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
 
-        const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
+        const Scalar velocity = elemFaceVars[scvf].velocitySelf();
 
         const bool insideIsUpstream = sign(scvf.outerNormalScalar()) == sign(velocity);
         const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars;
diff --git a/dumux/implicit/staggered/localresidual.hh b/dumux/implicit/staggered/localresidual.hh
index bdac723a8354a13d16fd7b9a6e83a628a821f85b..eca59a79aadd4fb6f1deba1e204c4e073f91003c 100644
--- a/dumux/implicit/staggered/localresidual.hh
+++ b/dumux/implicit/staggered/localresidual.hh
@@ -29,6 +29,11 @@
 
 namespace Dumux
 {
+
+namespace Properties
+{
+    NEW_PROP_TAG(ElementFaceVariables);
+}
 /*!
  * \ingroup CCModel
  * \ingroup StaggeredLocalResidual
@@ -58,7 +63,6 @@ class StaggeredLocalResidual
     using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
     using CellCenterSolutionVector = typename GET_PROP_TYPE(TypeTag, CellCenterSolutionVector);
     using FaceSolutionVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
-    using GlobalFaceVars = typename GET_PROP_TYPE(TypeTag, GlobalFaceVars);
     using CellCenterPrimaryVariables = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
     using FacePrimaryVariables = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
     using SolutionVector = typename GET_PROP_TYPE(TypeTag, SolutionVector);
@@ -66,6 +70,7 @@ class StaggeredLocalResidual
     using CellCenterResidual = typename GET_PROP_TYPE(TypeTag, CellCenterPrimaryVariables);
     using FaceResidual = typename GET_PROP_TYPE(TypeTag, FacePrimaryVariables);
     using FaceResidualVector = typename GET_PROP_TYPE(TypeTag, FaceSolutionVector);
+    using ElementFaceVariables = typename GET_PROP_TYPE(TypeTag, ElementFaceVariables);
 
 
     using DofTypeIndices = typename GET_PROP(TypeTag, DofTypeIndices);
@@ -118,8 +123,8 @@ public:
                         const FVElementGeometry& fvGeometry,
                         const ElementVolumeVariables& prevElemVolVars,
                         const ElementVolumeVariables& curElemVolVars,
-                        const GlobalFaceVars& prevFaceVars,
-                        const GlobalFaceVars& curFaceVars,
+                        const ElementFaceVariables& prevElemFaceVars,
+                        const ElementFaceVariables& curElemFaceVars,
                         const ElementBoundaryTypes &bcTypes,
                         const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
@@ -128,9 +133,9 @@ public:
         residual = 0.0;
         // ccStorageTerm_ = 0.0;
 
-        asImp_().evalVolumeTermForCellCenter_(residual, problem, element, fvGeometry, prevElemVolVars, curElemVolVars, prevFaceVars, curFaceVars, bcTypes);
-        asImp_().evalFluxesForCellCenter_(residual, problem, element, fvGeometry, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
-        asImp_().evalBoundaryForCellCenter_(residual, problem, element, fvGeometry, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
+        asImp_().evalVolumeTermForCellCenter_(residual, problem, element, fvGeometry, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, bcTypes);
+        asImp_().evalFluxesForCellCenter_(residual, problem, element, fvGeometry, curElemVolVars, curElemFaceVars, bcTypes, elemFluxVarsCache);
+        asImp_().evalBoundaryForCellCenter_(residual, problem, element, fvGeometry, curElemVolVars, curElemFaceVars, bcTypes, elemFluxVarsCache);
 
         return residual;
     }
@@ -157,17 +162,17 @@ public:
                   const SubControlVolumeFace& scvf,
                   const ElementVolumeVariables& prevElemVolVars,
                   const ElementVolumeVariables& curElemVolVars,
-                  const GlobalFaceVars& prevFaceVars,
-                  const GlobalFaceVars& curFaceVars,
+                  const ElementFaceVariables& prevElemFaceVars,
+                  const ElementFaceVariables& curElemFaceVars,
                   const ElementBoundaryTypes &bcTypes,
                   const ElementFluxVariablesCache& elemFluxVarsCache,
                   const bool resizeResidual = false) const
     {
         FaceResidual residual;
 
-        asImp_().evalVolumeTermForFace_(residual, problem, element, fvGeometry, scvf, prevElemVolVars, curElemVolVars, prevFaceVars, curFaceVars, bcTypes);
-        asImp_().evalFluxesForFace_(residual, problem, element, fvGeometry, scvf, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
-        asImp_().evalBoundaryForFace_(residual, problem, element, fvGeometry, scvf, curElemVolVars, curFaceVars, bcTypes, elemFluxVarsCache);
+        asImp_().evalVolumeTermForFace_(residual, problem, element, fvGeometry, scvf, prevElemVolVars, curElemVolVars, prevElemFaceVars, curElemFaceVars, bcTypes);
+        asImp_().evalFluxesForFace_(residual, problem, element, fvGeometry, scvf, curElemVolVars, curElemFaceVars, bcTypes, elemFluxVarsCache);
+        asImp_().evalBoundaryForFace_(residual, problem, element, fvGeometry, scvf, curElemVolVars, curElemFaceVars, bcTypes, elemFluxVarsCache);
 
         return residual;
     }
@@ -205,14 +210,14 @@ protected:
                                   const Element& element,
                                   const FVElementGeometry& fvGeometry,
                                   const ElementVolumeVariables& elemVolVars,
-                                  const GlobalFaceVars& faceVars,
+                                  const ElementFaceVariables& elemFaceVars,
                                   const ElementBoundaryTypes& bcTypes,
                                   const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
         for (auto&& scvf : scvfs(fvGeometry))
         {
             if(!scvf.boundary())
-                residual += asImp_().computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, faceVars, scvf, elemFluxVarsCache);
+                residual += asImp_().computeFluxForCellCenter(problem, element, fvGeometry, elemVolVars, elemFaceVars, scvf, elemFluxVarsCache);
         }
     }
 
@@ -225,12 +230,12 @@ protected:
                             const FVElementGeometry& fvGeometry,
                             const SubControlVolumeFace& scvf,
                             const ElementVolumeVariables& elemVolVars,
-                            const GlobalFaceVars& globalFaceVars,
+                            const ElementFaceVariables& elemFaceVars,
                             const ElementBoundaryTypes& bcTypes,
                             const ElementFluxVariablesCache& elemFluxVarsCache) const
     {
         if(!scvf.boundary())
-            residual += asImp_().computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, globalFaceVars, elemFluxVarsCache);
+            residual += asImp_().computeFluxForFace(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, elemFluxVarsCache);
     }
 
      /*!
@@ -240,7 +245,7 @@ protected:
                        const Element& element,
                        const FVElementGeometry& fvGeometry,
                        const ElementVolumeVariables& elemVolVars,
-                       const GlobalFaceVars& faceVars,
+                       const ElementFaceVariables& elemFaceVars,
                        const ElementBoundaryTypes& bcTypes,
                        const ElementFluxVariablesCache& elemFluxVarsCache)
     {
@@ -260,8 +265,8 @@ protected:
                                  const FVElementGeometry& fvGeometry,
                                  const ElementVolumeVariables& prevElemVolVars,
                                  const ElementVolumeVariables& curElemVolVars,
-                                 const GlobalFaceVars& prevFaceVars,
-                                 const GlobalFaceVars& curFaceVars,
+                                 const ElementFaceVariables& prevFaceVars,
+                                 const ElementFaceVariables& curFaceVars,
                                  const ElementBoundaryTypes &bcTypes) const
     {
         for(auto&& scv : scvs(fvGeometry))
@@ -287,8 +292,8 @@ protected:
                            const SubControlVolumeFace& scvf,
                            const ElementVolumeVariables& prevElemVolVars,
                            const ElementVolumeVariables& curElemVolVars,
-                           const GlobalFaceVars& prevFaceVars,
-                           const GlobalFaceVars& curFaceVars,
+                           const ElementFaceVariables& prevFaceVars,
+                           const ElementFaceVariables& curFaceVars,
                            const ElementBoundaryTypes &bcTypes) const
     {
         // the source term:
@@ -310,8 +315,8 @@ protected:
                                  const FVElementGeometry& fvGeometry,
                                  const ElementVolumeVariables& prevElemVolVars,
                                  const ElementVolumeVariables& curElemVolVars,
-                                 const GlobalFaceVars& prevFaceVars,
-                                 const GlobalFaceVars& curFaceVars,
+                                 const ElementFaceVariables& prevFaceVars,
+                                 const ElementFaceVariables& curFaceVars,
                                  const ElementBoundaryTypes &bcTypes) const
     {
         for(auto&& scv : scvs(fvGeometry))
@@ -361,8 +366,8 @@ protected:
                            const SubControlVolumeFace& scvf,
                            const ElementVolumeVariables& prevElemVolVars,
                            const ElementVolumeVariables& curElemVolVars,
-                           const GlobalFaceVars& prevFaceVars,
-                           const GlobalFaceVars& curFaceVars,
+                           const ElementFaceVariables& prevFaceVars,
+                           const ElementFaceVariables& curFaceVars,
                            const ElementBoundaryTypes &bcTypes) const
     {
         const auto& scv = fvGeometry.scv(scvf.insideScvIdx());