diff --git a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
index bbd5e1981277bc3b39973cbc79608c72b1fa8246..d1410815f2503d731f371af430d74a72a09250d6 100644
--- a/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/fluxvariables.hh
@@ -547,7 +547,8 @@ private:
                                                                     const FaceVariables& faceVars,
                                                                     const std::optional<BoundaryTypes>& currentScvfBoundaryTypes,
                                                                     const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes,
-                                                                    const int localSubFaceIdx)
+                                                                    const int localSubFaceIdx,
+                                                                    bool useZeroNormalGradient = false)
     {
         const auto eIdx = scvf.insideScvIdx();
         const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
@@ -568,8 +569,12 @@ private:
                              ? insideVolVars.effectiveViscosity()
                              : (insideVolVars.effectiveViscosity() + outsideVolVars.effectiveViscosity()) * 0.5;
 
+
+        if(currentScvfBoundaryTypes->isDirichlet(Indices::pressureIdx))
+            useZeroNormalGradient = true;
+
         // Consider the shear stress caused by the gradient of the velocities normal to our face of interest.
-        if (!enableUnsymmetrizedVelocityGradient)
+        if (!enableUnsymmetrizedVelocityGradient && !useZeroNormalGradient)
         {
             if (!scvf.boundary() ||
                 currentScvfBoundaryTypes->isDirichlet(Indices::velocity(lateralFace.directionIndex())) ||
@@ -584,11 +589,8 @@ private:
         // Consider the shear stress caused by the gradient of the velocities parallel to our face of interest.
         // If we have a Dirichlet condition for the pressure at the lateral face we assume to have a zero velocityGrad_ij velocity gradient
         // so we can skip the computation.
-        if (!lateralFace.boundary() || !lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx))
-        {
-            const Scalar velocityGrad_ij = VelocityGradients::velocityGradIJ(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
-            lateralDiffusiveFlux -= muAvg * velocityGrad_ij * lateralFace.directionSign();
-        }
+        const Scalar velocityGrad_ij = VelocityGradients::velocityGradIJ(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
+        lateralDiffusiveFlux -= muAvg * velocityGrad_ij * lateralFace.directionSign();
 
         // Account for the area of the staggered lateral face (0.5 of the coinciding scfv).
         return lateralDiffusiveFlux * lateralFace.area() * 0.5 * extrusionFactor_(elemVolVars, lateralFace);
diff --git a/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh b/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh
index 8993245d319c321dfff1db59811609ae23d86475..e8232cee3b1c78f48c374754085a4012fdec10f2 100644
--- a/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh
+++ b/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh
@@ -521,16 +521,13 @@ private:
                                                    const int localSubFaceIdx)
     {
         // Find out what boundary type is set on the lateral face
-        const bool useZeroGradient = lateralFaceBoundaryTypes && (lateralFaceBoundaryTypes->isSymmetry()
-                                                                  || lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx));
         const bool lateralFaceHasBJS = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex()));
         const bool lateralFaceHasDirichletVelocity = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex()));
-        const Scalar velocitySelf = faceVars.velocitySelf();
 
         // If there is a Dirichlet condition for the pressure we assume zero gradient for the velocity,
         // so the velocity at the boundary equal to that on the scvf.
-        if (useZeroGradient)
-            return velocitySelf;
+        if (!lateralFaceHasBJS && !lateralFaceHasDirichletVelocity)
+            return faceVars.velocitySelf();
 
         if (lateralFaceHasBJS)
             return VelocityGradients::beaversJosephVelocityAtLateralScvf(problem, element, fvGeometry, scvf,  faceVars,
diff --git a/dumux/freeflow/navierstokes/staggered/velocitygradients.hh b/dumux/freeflow/navierstokes/staggered/velocitygradients.hh
index a9c3c452b31e9d5ed9ba383d29cf3a592bd87986..97dfe5fa744998f1727d60a22ebe5caa50ac55c4 100644
--- a/dumux/freeflow/navierstokes/staggered/velocitygradients.hh
+++ b/dumux/freeflow/navierstokes/staggered/velocitygradients.hh
@@ -127,7 +127,7 @@ public:
                                                           currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
             }
             else
-                DUNE_THROW(Dune::InvalidStateException, "Invalid lateral boundary type at " << lateralScvf.center());
+                return innerParallelVelocity;
         }();
 
         // The velocity gradient already accounts for the orientation
@@ -181,10 +181,6 @@ public:
         const auto eIdx = scvf.insideScvIdx();
         const auto& lateralScvf = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx);
 
-        // Assume a zero velocity gradient for pressure boundary conditions.
-        if (currentScvfBoundaryTypes && currentScvfBoundaryTypes->isDirichlet(Indices::pressureIdx))
-            return 0.0;
-
         // For the velocityGrad_ji gradient, get the velocities perpendicular to the velocity at the current scvf.
         // The inner one is located at staggered face within the own element,
         // the outer one at the respective staggered face of the element on the other side of the
@@ -206,7 +202,7 @@ public:
                                                           currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
             }
             else
-                DUNE_THROW(Dune::InvalidStateException, "Invalid lateral boundary types at " << lateralScvf.center());
+                return innerLateralVelocity;
         }();
 
         // Calculate the velocity gradient in positive coordinate direction.
@@ -261,15 +257,10 @@ public:
 
             if (unsymmetrizedGradientForBJ)
                 return 0.0;
-
-            if (lateralScvf.boundary())
-            {
-                if (lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx) ||
-                    lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex())))
-                    return 0.0;
-            }
-
-            return velocityGradIJ(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
+            else if (lateralScvf.boundary() && lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex())))
+                return 0.0;
+            else
+                return velocityGradIJ(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
         }();
 
         return problem.beaversJosephVelocity(element,
@@ -321,15 +312,10 @@ public:
 
             if (unsymmetrizedGradientForBJ)
                 return 0.0;
-
-            if (scvf.boundary())
-            {
-                if (currentScvfBoundaryTypes->isDirichlet(Indices::pressureIdx) ||
-                    currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex())))
-                    return 0.0;
-            }
-
-            return velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
+            else if (scvf.boundary() && currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralScvf.directionIndex())))
+                return 0.0;
+            else
+                return velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx);
         }();
 
         return problem.beaversJosephVelocity(element,