diff --git a/dumux/discretization/staggered/freeflow/fickslaw.hh b/dumux/discretization/staggered/freeflow/fickslaw.hh
index dd9963ded86cb3cdde99f08666075286efec4785..cd70d075d99811bee364e9e7149524fe17431d02 100644
--- a/dumux/discretization/staggered/freeflow/fickslaw.hh
+++ b/dumux/discretization/staggered/freeflow/fickslaw.hh
@@ -100,7 +100,7 @@ public:
 
         const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
         const auto& insideVolVars = elemVolVars[insideScv];
-        const auto& outsideVolVars = scvf.boundary() ?  insideVolVars : elemVolVars[scvf.outsideScvIdx()];
+        const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
 
         const Scalar insideMolarDensity = insideVolVars.molarDensity();
 
@@ -117,27 +117,14 @@ public:
             if(scvf.boundary())
             {
                 const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
-                if(bcTypes.isOutflow(eqIdx))
+                if(bcTypes.isOutflow(momentumBalanceIdx) || bcTypes.isNeumann(eqIdx)) // TODO: catch neumann and outflow in localResidual's evalBoundary_()
                     return flux;
-                else if(bcTypes.isNeumann(eqIdx))
-                    return flux; // TODO: implement neumann
-                else
-                {
-                    // get the Dirichlet value for the mole fraction on the boundary
-                    Scalar dirichletMoleFraction = problem.dirichletAtPos(scvf.center())[eqIdx];
-                    // convert value to a mole fraction if user specifies a mass fraction
-                    if(!useMoles)
-                        dirichletMoleFraction /= FluidSystem::molarMass(compIdx) / insideVolVars.fluidState().averageMolarMass(phaseIdx);
-                    flux[compIdx] = insideMolarDensity * tij * (insideMoleFraction - dirichletMoleFraction);
-                }
-            }
-            else
-            {
-                const Scalar outsideMolarDensity = outsideVolVars.molarDensity();
-                const Scalar avgDensity = 0.5*(insideMolarDensity + outsideMolarDensity);
-                const Scalar outsideMoleFraction = outsideVolVars.moleFraction(phaseIdx, compIdx);
-                flux[compIdx] = avgDensity * tij * (insideMoleFraction - outsideMoleFraction);
             }
+
+            const Scalar outsideMolarDensity = scvf.boundary() ? insideVolVars.molarDensity() : outsideVolVars.molarDensity();
+            const Scalar avgDensity = 0.5*(insideMolarDensity + outsideMolarDensity);
+            const Scalar outsideMoleFraction = outsideVolVars.moleFraction(phaseIdx, compIdx);
+            flux[compIdx] = avgDensity * tij * (insideMoleFraction - outsideMoleFraction);
             ++eqIdx;
         }
 
@@ -179,9 +166,7 @@ public:
     {
         Scalar tij = 0.0;
         const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx());
-        const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
         const auto& insideVolVars = elemVolVars[insideScv];
-        const auto& outsideVolVars = scvf.boundary() ?  insideVolVars : elemVolVars[scvf.outsideScvIdx()];
 
         const Scalar insideDistance = (insideScv.dofPosition() - scvf.ipGlobal()).two_norm();
         const Scalar insideD = insideVolVars.diffusionCoefficient(phaseIdx, compIdx);
@@ -191,6 +176,8 @@ public:
             tij = scvf.area() * ti;
         else
         {
+            const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx());
+            const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
             const Scalar outsideDistance = (outsideScv.dofPosition() - scvf.ipGlobal()).two_norm();
             const Scalar outsideD = outsideVolVars.diffusionCoefficient(phaseIdx, compIdx);
             const Scalar tj = calculateOmega_(outsideDistance, outsideD, 1.0);
diff --git a/dumux/freeflow/staggerednc/fluxvariables.hh b/dumux/freeflow/staggerednc/fluxvariables.hh
index da1f309a1b1d1bc24fc78662ab1599d4f66ddccb..0083a00dcbf0f06274187e1308751bd28decbbcd 100644
--- a/dumux/freeflow/staggerednc/fluxvariables.hh
+++ b/dumux/freeflow/staggerednc/fluxvariables.hh
@@ -127,7 +127,16 @@ private:
         const auto& insideVolVars = elemVolVars[insideScv];
 
         // if we are on an inflow/outflow boundary, use the volVars of the element itself
-        const auto& outsideVolVars = scvf.boundary() ?  insideVolVars : elemVolVars[scvf.outsideScvIdx()];
+        // TODO: catch neumann and outflow in localResidual's evalBoundary_()
+        bool isOutflow = false;
+        if(scvf.boundary())
+        {
+            const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
+                if(bcTypes.isOutflow(momentumBalanceIdx))
+                    isOutflow = true;
+        }
+
+        const auto& outsideVolVars = isOutflow ? insideVolVars : elemVolVars[scvf.outsideScvIdx()];
 
         const Scalar velocity = globalFaceVars.faceVars(scvf.dofIndex()).velocity();
 
@@ -142,34 +151,20 @@ private:
         for (int compIdx = 0; compIdx < numComponents; ++compIdx)
         {
             // get equation index
-            auto eqIdx = conti0EqIdx + compIdx;
+            const auto eqIdx = conti0EqIdx + compIdx;
+            if (eqIdx == replaceCompEqIdx)
+                continue;
 
-            const Scalar upstreamDensity = useMoles ? upstreamVolVars.molarDensity() : upstreamVolVars.density();
-            const Scalar downstreamDensity = useMoles ? downstreamVolVars.molarDensity() : downstreamVolVars.density();
             const Scalar upstreamFraction = useMoles ? upstreamVolVars.moleFraction(phaseIdx, compIdx) : upstreamVolVars.massFraction(phaseIdx, compIdx);
             const Scalar downstreamFraction = useMoles ? downstreamVolVars.moleFraction(phaseIdx, compIdx) : downstreamVolVars.massFraction(phaseIdx, compIdx);
 
-            Scalar advFlux = 0.0;
-
-            if(scvf.boundary() && eqIdx > conti0EqIdx)
-            {
-                const auto bcTypes = problem.boundaryTypesAtPos(scvf.center());
-                if(bcTypes.isDirichlet(eqIdx))
-                    advFlux = upstreamDensity * problem.dirichletAtPos(scvf.center())[eqIdx] * velocity;
-                if(bcTypes.isOutflow(eqIdx))
-                    advFlux = upstreamDensity * upstreamFraction * velocity;
-            }
-            else
-                advFlux = (upWindWeight * upstreamDensity * upstreamFraction +
-                          (1.0 - upWindWeight) * downstreamDensity * downstreamFraction) * velocity;
-
-            if (eqIdx != replaceCompEqIdx)
-                flux[eqIdx] += advFlux;
-
-            // in case one balance is substituted by the total mass balance
-            if (replaceCompEqIdx < numComponents)
-                flux[replaceCompEqIdx] += advFlux;
+            flux[eqIdx] = (upWindWeight * upstreamDensity * upstreamFraction +
+                          (1.0 - upWindWeight) * downstreamDensity * downstreamFraction)
+                          * velocity;
         }
+        // in case one balance is substituted by the total mass balance
+        if (replaceCompEqIdx < numComponents)
+            flux[replaceCompEqIdx] = (upWindWeight * upstreamDensity + (1.0 - upWindWeight) * downstreamDensity) * velocity;
 
         flux *= scvf.area() * sign(scvf.outerNormalScalar());
         return flux;
diff --git a/dumux/freeflow/staggerednc/localresidual.hh b/dumux/freeflow/staggerednc/localresidual.hh
index b96ff69a914d4d130bd84fe44b19fb0eca2da1f3..0a8786a7c96cef4b2dcc3490ecf8c8dbbde9dc87 100644
--- a/dumux/freeflow/staggerednc/localresidual.hh
+++ b/dumux/freeflow/staggerednc/localresidual.hh
@@ -119,7 +119,7 @@ class StaggeredNavierStokesResidualImpl<TypeTag, true, false> : public Staggered
         }
             // in case one balance is substituted by the total mass balance
             if(replaceCompEqIdx < numComponents)
-                storage[replaceCompEqIdx] += density;
+                storage[replaceCompEqIdx] = density;
 
         //TODO: energy balance
         return storage;