diff --git a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
index a54d0e24acb910cd95f9d9ef5ff9351b1997e68e..1022cf31eda012dd313a34bb8294300513b0289c 100644
--- a/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
+++ b/dumux/discretization/cellcentered/mpfa/fluxvariablescachefiller.hh
@@ -58,6 +58,7 @@ class CCMpfaFluxVariablesCacheFillerImplementation<TypeTag, true, false, false>
     using Element = typename GridView::template Codim<0>::Entity;
 
     static const int numEq = GET_PROP_VALUE(TypeTag, NumEq);
+    static const bool useTpfaBoundary = GET_PROP_VALUE(TypeTag, UseTpfaBoundary);
     static const bool solDependentParams = GET_PROP_VALUE(TypeTag, SolutionDependentParameters);
 
 public:
@@ -68,36 +69,46 @@ public:
                                  const FVElementGeometry& fvGeometry,
                                  const ElementVolumeVariables& elemVolVars,
                                  const SubControlVolumeFace& scvf,
-                                 FluxVarsCacheVector& fluxVarsCache)
+                                 FluxVarsCacheVector& fluxVarsCache,
+                                 const bool updateNeumannOnly = false)
     {
+        // if we only want to update the neumann fluxes, skip the rest if we don't touch the boundary
+        const bool boundary = problem.model().globalFvGeometry().scvfTouchesBoundary(scvf);
+        if (updateNeumannOnly && !boundary)
+            return;
+
         // lambda function to get the permeability tensor
         const auto* prob = &problem;
-        auto permFunction = [prob](const Element& element, const VolumeVariables& volVars, const SubControlVolume& scv)
+        auto permFunction = [prob](const Element& element,
+                                   const VolumeVariables& volVars,
+                                   const SubControlVolume& scv)
                             { return prob->spatialParams().intrinsicPermeability(scv, volVars); };
 
-        // update the flux var caches for this scvf
-        if (problem.model().globalFvGeometry().scvfTouchesBoundary(scvf))
+        // get the interaction volume seed and update flux var caches
+        const auto& seed = boundary ? problem.model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf) : problem.model().globalFvGeometry().interactionVolumeSeed(scvf);
+        if (boundary)
         {
-            const auto& boundarySeed = problem.model().globalFvGeometry().boundaryInteractionVolumeSeed(scvf);
-            BoundaryInteractionVolume iv(boundarySeed, problem, fvGeometry, elemVolVars);
+            BoundaryInteractionVolume iv(seed, problem, fvGeometry, elemVolVars);
             iv.solveLocalSystem(permFunction);
 
             // we assume phaseIdx = eqIdx here for purely advective problems
             for (unsigned int eqIdx = 0; eqIdx < numEq; ++eqIdx)
             {
                 // lambda function defining the upwind factor of the advective flux
-                auto advectionUpwindFunction = [eqIdx](const VolumeVariables& volVars) { return volVars.density(eqIdx)/volVars.viscosity(eqIdx); };
+                auto advectionUpwindFunction = [eqIdx](const VolumeVariables& volVars)
+                                               { return volVars.density(eqIdx)/volVars.viscosity(eqIdx); };
                 iv.assembleNeumannFluxes(advectionUpwindFunction, eqIdx);
 
-                // update flux variables cache, only the neumann fluxes have to be set for each phase
+                // update flux variables cache for the scvf
                 const auto scvfIdx = scvf.index();
                 auto& cache = fluxVarsCache[scvfIdx];
-                if (eqIdx == 0)
+                cache.updatePhaseNeumannFlux(problem, element, fvGeometry, elemVolVars, scvf, iv, eqIdx);
+                if (eqIdx == numEq-1)
                 {
-                    cache.updateAdvection(problem, element, fvGeometry, elemVolVars, scvf, iv);
+                    if (!updateNeumannOnly)
+                        cache.updateAdvection(problem, element, fvGeometry, elemVolVars, scvf, iv);
                     cache.setUpdated();
                 }
-                cache.updatePhaseNeumannFlux(problem, element, fvGeometry, elemVolVars, scvf, iv, eqIdx);
 
                 // update flux variable caches of the other scvfs of the interaction volume
                 for (const auto scvfIdxJ : iv.globalScvfs())
@@ -107,19 +118,19 @@ public:
                         const auto& scvfJ = fvGeometry.scvf(scvfIdxJ);
                         const auto elementJ = problem.model().globalFvGeometry().element(scvfJ.insideScvIdx());
                         auto& cacheJ = fluxVarsCache[scvfIdxJ];
-                        if (eqIdx == 0)
+                        cacheJ.updatePhaseNeumannFlux(problem, elementJ, fvGeometry, elemVolVars, scvfJ, iv, eqIdx);
+                        if (eqIdx == numEq-1)
                         {
-                            cacheJ.updateAdvection(problem, elementJ, fvGeometry, elemVolVars, scvfJ, iv);
+                            if (!updateNeumannOnly)
+                                cacheJ.updateAdvection(problem, elementJ, fvGeometry, elemVolVars, scvfJ, iv);
                             cacheJ.setUpdated();
                         }
-                        cacheJ.updatePhaseNeumannFlux(problem, elementJ, fvGeometry, elemVolVars, scvfJ, iv, eqIdx);
                     }
                 }
             }
         }
         else
         {
-            const auto& seed = problem.model().globalFvGeometry().interactionVolumeSeed(scvf);
             InteractionVolume iv(seed, problem, fvGeometry, elemVolVars);
             iv.solveLocalSystem(permFunction);
 
@@ -154,7 +165,13 @@ public:
                                    FluxVarsCacheVector& fluxVarsCache)
     {
         if (!solDependentParams)
+        {
             fluxVarsCache[scvf.index()].setUpdated();
+
+            // if we do not use tpfa on the boundary, we have to update the neumann fluxes
+            if (!useTpfaBoundary)
+                fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCache, true);
+        }
         else
             fillFluxVarCache(problem, element, fvGeometry, elemVolVars, scvf, fluxVarsCache);
     }
diff --git a/dumux/porousmediumflow/implicit/fluxvariablescache.hh b/dumux/porousmediumflow/implicit/fluxvariablescache.hh
index b82c3191b6d78176644287b3d774ca010eb0ed42..3e7de05f2ea6703e84a35d8fb17d97a6850854d5 100644
--- a/dumux/porousmediumflow/implicit/fluxvariablescache.hh
+++ b/dumux/porousmediumflow/implicit/fluxvariablescache.hh
@@ -175,7 +175,6 @@ class PorousMediumMpfaFluxVariablesCache<TypeTag, true, false, false>
     using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
     using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
     using BoundaryInteractionVolume = typename GET_PROP_TYPE(TypeTag, BoundaryInteractionVolume);
-    using InteractionVolume = typename GET_PROP_TYPE(TypeTag, InteractionVolume);
     using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
     using AdvectionType = typename GET_PROP_TYPE(TypeTag, AdvectionType);
     using Element = typename GridView::template Codim<0>::Entity;
@@ -195,11 +194,13 @@ public:
     // the constructor
     PorousMediumMpfaFluxVariablesCache() : isUpdated_(false)
     {
+        // We have to initialize to zero (for inner interaction volumes)
         for (auto& nFlux : phaseNeumannFluxes_)
             nFlux = 0.0;
     }
 
     // update cached objects
+    template<typename InteractionVolume>
     void updateAdvection(const Problem& problem,
                          const Element& element,
                          const FVElementGeometry& fvGeometry,
@@ -222,6 +223,7 @@ public:
         tij_ = interactionVolume.getTransmissibilities(localIndexPair);
     }
 
+    template<typename InteractionVolume>
     void updatePhaseNeumannFlux(const Problem& problem,
                                 const Element& element,
                                 const FVElementGeometry& fvGeometry,