diff --git a/dune/dorie/solver/operator_DG.hh b/dune/dorie/solver/operator_DG.hh
index 9253f9aa0bec6be512be4406c035ae9f4c400d96..d4708381e4af942df168fbf76f19bb90d0920353 100755
--- a/dune/dorie/solver/operator_DG.hh
+++ b/dune/dorie/solver/operator_DG.hh
@@ -474,10 +474,22 @@ public:
           // compute numerical flux estimation
           const RF numFlux = satCond_s * ( - (gradu_s * normal) + penalty * jump);
 
+          RF head_upwind;
+          if (numFlux > 0)
+            head_upwind = param.headMillerToRef(u_s, p_global_s);
+          else
+            head_upwind = g; // Miller here?!
+
+          // compute saturation from matrix head
+          const RF saturation_upwind = param.saturation(head_upwind, p_global_s);
+
+          // compute relative conductivity
+          const RF relCond_upwind = param.condFactor(saturation_upwind, p_global_s);
+
           // switches between Neumann and Dirichlet bc
           // this choice is kept constant for one time step
           normal_flux = boundary.j(ig.intersection(),it.position(),time);
-          if(numFlux < normal_flux) {
+          if(relCond_upwind * numFlux < normal_flux) {
             bcType = BCType::dirichlet;
             bc_type_cache.emplace(it_id,BoundaryCondition::Dirichlet);
           }
@@ -557,7 +569,7 @@ public:
           // consistency term
           // + penalty term
           for (unsigned int i = 0; i<lfsv_s.size(); i++)
-            r_s.accumulate(lfsv_s,i, numFlux * phi_s[i] * factor);
+            r_s.accumulate(lfsv_s,i, relCond_upwind * numFlux * phi_s[i] * factor);
 
           // update residual (symmetry term)
           // (non-)symmetric IP term