diff --git a/dune/dorie/solver/operator_DG.hh b/dune/dorie/solver/operator_DG.hh
index 6e71fd75b8b932a06670559157dd68e57a45d30b..8b26dd56e789a1671cf9ee0c0184b71e6a13b26a 100755
--- a/dune/dorie/solver/operator_DG.hh
+++ b/dune/dorie/solver/operator_DG.hh
@@ -57,7 +57,7 @@ template<typename Traits, typename Parameter, typename Boundary, typename Source
   :
   // public Dune::PDELab::NumericalJacobianVolume<RichardsDGSpatialOperator<Traits, Parameter, Boundary, SourceTerm, FEM, adjoint> >,
   // public Dune::PDELab::NumericalJacobianSkeleton<RichardsDGSpatialOperator<Traits, Parameter, Boundary, SourceTerm, FEM, adjoint> >,
-  public Dune::PDELab::NumericalJacobianBoundary<RichardsDGSpatialOperator<Traits, Parameter, Boundary, SourceTerm, FEM, adjoint> >,
+  // public Dune::PDELab::NumericalJacobianBoundary<RichardsDGSpatialOperator<Traits, Parameter, Boundary, SourceTerm, FEM, adjoint> >,
   public Dune::PDELab::FullSkeletonPattern,
   public Dune::PDELab::FullVolumePattern,
   public Dune::PDELab::LocalOperatorDefaultFlags,
@@ -130,7 +130,7 @@ public:
     :
     // Dune::PDELab::NumericalJacobianVolume<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
     // Dune::PDELab::NumericalJacobianSkeleton<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
-    Dune::PDELab::NumericalJacobianBoundary<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
+    // Dune::PDELab::NumericalJacobianBoundary<RichardsDGSpatialOperator<Traits,Parameter,Boundary,SourceTerm,FEM,adjoint> >(1.e-7),
     param(param_), boundary(boundary_), sourceTerm(sourceTerm_), method(method_), weights(weights_),
     penalty_factor(config.get<RF>("dg.penaltyFactor")), mapper(view_),
     intorderadd(intorderadd_), quadrature_factor(quadrature_factor_),
@@ -576,14 +576,14 @@ public:
           mat_ss.accumulate(lfsv_s,i,lfsu_s,j,
             relCond_upwind 
             * ( omega_s * satCond_s
-                * ( - ( tgradphi_s[j] * normal ) + penalty * phi_s[j] ) 
+                * ( - ( tgradphi_s[j] * normal ) + penalty * phi_s[j] )
               + omega_n * satCond_n * penalty * phi_s[j] )
             * phi_s[i] * factor
           );
           // differentiate symmetry term
           mat_ss.accumulate(lfsv_s,i,lfsu_s,j,
-            theta * omega_s * (tgradphi_s[i] * normal) * satCond_s
-            * phi_s[j] * ( rel_cond_upwind_deriv_s + 1 )
+            theta * omega_s * satCond_s * (tgradphi_s[i] * normal)
+            * phi_s[j] * ( rel_cond_upwind_deriv_s + 1.0 )
             * factor
           );
         }
@@ -604,7 +604,7 @@ public:
           // differentiate symmetry term
           mat_nn.accumulate(lfsv_n,i,lfsu_n,j,
             theta * omega_n * (tgradphi_n[i] * normal) * satCond_n
-            * phi_n[j] * ( rel_cond_upwind_deriv_n - 1 )
+            * phi_n[j] * ( rel_cond_upwind_deriv_n - 1.0 )
             * factor
           );
         }
@@ -625,7 +625,7 @@ public:
           // differentiate symmetry term
           mat_sn.accumulate(lfsv_s,i,lfsu_n,j,
             theta * omega_s * (tgradphi_s[i] * normal) * satCond_s
-            * phi_n[j] * ( rel_cond_upwind_deriv_n - 1 )
+            * phi_n[j] * ( rel_cond_upwind_deriv_n - 1.0 )
             * factor
           );
         }
@@ -646,7 +646,7 @@ public:
           // differentiate symmetry term
           mat_ns.accumulate(lfsv_n,i,lfsu_s,j,
             theta * omega_n * (tgradphi_n[i] * normal) * satCond_n
-            * phi_s[j] * ( rel_cond_upwind_deriv_s + 1 )
+            * phi_s[j] * ( rel_cond_upwind_deriv_s + 1.0 )
             * factor
           );
         }
@@ -899,6 +899,223 @@ public:
       }
     }
 
+  template<typename IG, typename LFSU, typename X, typename LFSV, typename LocalMatrix>
+  void jacobian_boundary (const IG& ig,
+    const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s,
+    LocalMatrix& mat_ss) const
+  {
+    // get polynomial degree
+    const int order_s  = lfsu_s.finiteElement().localBasis().order();
+    const int degree   = order_s;
+    const int intorder = intorderadd + quadrature_factor * degree;
+
+    // geometric factor of penalty
+    const RF h_F = ig.inside().geometry().volume()/ig.geometry().volume(); // Houston!
+
+    // get element geometry
+    auto gtface = ig.geometry();
+
+    // loop over quadrature points and integrate normal flux
+    std::size_t i = 0;
+    for (const auto& it : quadratureRule(gtface,intorder))
+    {
+      // calculate unique id
+      const auto it_id = create_unique_id(ig.intersection(),i);
+      i++;
+
+      // position of quadrature point in local coordinates of elements
+      const auto p_local_s = ig.geometryInInside().global(it.position());
+
+      // evaluate basis functions
+      const auto& phi_s = cache[order_s].evaluateFunction(p_local_s,lfsu_s.finiteElement().localBasis());
+
+      // integration factor
+      const RF factor = it.weight() * ig.geometry().integrationElement(it.position());
+
+      // evaluate u
+      RF u_s = 0.;
+      for (unsigned int i = 0; i<lfsu_s.size(); i++)
+        u_s += x_s(lfsu_s,i) * phi_s[i];
+
+      // check if a BC type is cached
+      BoundaryCondition::Type bc;
+      const auto bc_cached = bc_type_cache.find(it_id);
+      if(bc_cached!=bc_type_cache.end())
+        bc = bc_cached->second;
+      else // query type of boundary condition
+        bc = boundary.bc(ig.intersection(),it.position(),time);
+
+      // set the internal bcType flag
+      auto bcType = BCType::none;
+      RF normal_flux;
+
+      if (BoundaryCondition::isNeumann(bc))
+      {
+        // evaluate flux boundary condition
+        normal_flux = boundary.j(ig.intersection(),it.position(),time);
+
+        bcType = BCType::vonNeumann;
+      }
+
+      else if (BoundaryCondition::isDirichlet(bc))
+      {
+        bcType = BCType::dirichlet;
+      }
+
+      else if (BoundaryCondition::isLimitedInFlux(bc))
+      {
+        // modified Neumann boundary condition, set flux to 0 where potential is positive
+        normal_flux = .0;
+        if(u_s < .0) normal_flux = boundary.j(ig.intersection(),it.position(),time);
+
+        bcType = BCType::vonNeumann;
+      }
+
+      else if (BoundaryCondition::isEvaporation(bc))
+      {
+        // estimate the dirichlet flux
+        // evaluate gradient of basis functions (we assume Galerkin method lfsu = lfsv)
+        const auto& gradphi_s = cache[order_s].evaluateJacobian(p_local_s,lfsu_s.finiteElement().localBasis());
+
+        // transform gradients to real element
+        const auto jac = ig.inside().geometry().jacobianInverseTransposed(p_local_s);
+        std::vector<Vector> tgradphi_s(lfsu_s.size());
+        for (unsigned int i = 0; i<lfsu_s.size(); i++)
+          jac.mv(gradphi_s[i][0],tgradphi_s[i]);
+
+        // compute gradient of u
+        Vector gradu_s(0.);
+        for (unsigned int i = 0; i<lfsu_s.size(); i++)
+          gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
+
+        // update with gravity vector
+        gradu_s -= param.gravity();
+
+        // face normal vector
+        const Domain normal = ig.unitOuterNormal(it.position());
+
+        // jump relative to Dirichlet value
+        const RF g = boundary.g(ig.intersection(),it.position(),time);
+        const RF jump = u_s - g;
+
+        const Domain p_global_s = ig.inside().geometry().global(p_local_s);
+
+        // conductivity at saturation
+        const RF satCond_s = param.condRefToMiller(param.K(p_global_s), p_global_s);
+
+        const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);
+
+        // compute numerical flux estimation
+        const RF numFlux = satCond_s * ( - (gradu_s * normal) + penalty * jump);
+
+        // 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) {
+          bcType = BCType::dirichlet;
+          bc_type_cache.emplace(it_id,BoundaryCondition::Dirichlet);
+        }
+        else {
+          bcType = BCType::vonNeumann;
+          bc_type_cache.emplace(it_id,BoundaryCondition::Neumann);
+        }
+      }
+
+      else if (BoundaryCondition::isOther(bc))
+        bcType = BCType::none;
+
+
+      // update residual according to byType flag
+      if (bcType == BCType::vonNeumann)
+      {
+        // update jacobian (flux term)
+        for (unsigned int i = 0; i<lfsv_s.size(); ++i) {
+          for (unsigned int j = 0; j<lfsu_s.size(); ++j) {
+            mat_ss.accumulate(lfsv_s,i,lfsu_s,j, 0.0);
+          }
+        }
+
+        continue;
+      }
+
+      else if (bcType == BCType::dirichlet)
+      {
+        // evaluate gradient of basis functions (we assume Galerkin method lfsu = lfsv)
+        const auto& gradphi_s = cache[order_s].evaluateJacobian(p_local_s,lfsu_s.finiteElement().localBasis());
+
+        // transform gradients to real element
+        const auto jac = ig.inside().geometry().jacobianInverseTransposed(p_local_s);
+        std::vector<Vector> tgradphi_s(lfsu_s.size());
+        for (unsigned int i = 0; i<lfsu_s.size(); i++)
+          jac.mv(gradphi_s[i][0],tgradphi_s[i]);
+
+        // compute gradient of u
+        Vector gradu_s(0.);
+        for (unsigned int i = 0; i<lfsu_s.size(); i++)
+          gradu_s.axpy(x_s(lfsu_s,i),tgradphi_s[i]);
+
+        // update with gravity vector
+        gradu_s -= param.gravity();
+
+        // face normal vector
+        const Domain normal = ig.unitOuterNormal(it.position());
+
+        // jump relative to Dirichlet value
+        const RF g = boundary.g(ig.intersection(),it.position(),time);
+        const RF jump = u_s - g;
+
+        const Domain p_global_s = ig.inside().geometry().global(p_local_s);
+
+        // conductivity at saturation
+        const RF satCond_s = param.condRefToMiller(param.K(p_global_s), p_global_s);
+
+        const RF penalty = penalty_factor / h_F * degree * (degree + dim - 1);
+
+        // compute numerical flux
+        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);
+
+        // derivatives of relCond to s
+        RF rel_cond_upwind_deriv_s;
+        if (numFlux > 0){
+          // relCond depends on s
+          rel_cond_upwind_deriv_s = param.cond_factor_deriv(saturation_upwind, p_global_s);
+          rel_cond_upwind_deriv_s *= param.saturation_deriv(head_upwind,p_global_s);
+        }
+        else {
+          // relCond depends on boundary condition
+          rel_cond_upwind_deriv_s = 0.0;
+        }
+
+        for (unsigned int i = 0; i<lfsv_s.size(); ++i) {
+          for (unsigned int j = 0; j<lfsu_s.size(); ++j) {
+            // update jacobian (flux term)
+            mat_ss.accumulate(lfsv_s,i,lfsu_s,j, satCond_s * ( - ( tgradphi_s[j] * normal ) + penalty * phi_s[j] ) * phi_s[i] * factor );
+            // update jacobian (symmetry term)
+            mat_ss.accumulate(lfsv_s,i,lfsu_s,j, theta * ( tgradphi_s[i] * normal ) * satCond_s * phi_s[j] * ( rel_cond_upwind_deriv_s + 1.0 ) * factor );
+          }
+        }
+
+        continue;
+      }
+
+      else if (bcType == BCType::none)
+        continue;
+
+    } // quadrature rule
+  } // jacobian_boundary
+
   /**
    * @brief Volume integral depending only on test functions
    */