diff --git a/dune/dorie/model/transport/local_operator_FV.hh b/dune/dorie/model/transport/local_operator_FV.hh
index c36fa6fd430f87af2d138478cfe5e5ff88bae12f..5b8f411b69a4afd8bb5afe0897a0ea5912760ffd 100644
--- a/dune/dorie/model/transport/local_operator_FV.hh
+++ b/dune/dorie/model/transport/local_operator_FV.hh
@@ -124,8 +124,8 @@ public:
     , _boundary(boundary)
     , _dirichlet_mode(dirichlet_mode)
   {
-    assert(_gf_water_flux);
-    assert(_gf_water_content);
+    assert(_param);
+    assert(_boundary);
   }
 
   /*---------------------------------------------------------------------*//**
@@ -155,6 +155,8 @@ public:
     const LFSU& lfsu_o, const X& x_o, const LFSV& lfsv_o, 
     R& r_i, R& r_o) const
   {
+    assert_pointers();
+
     // Get entity dim
     const int dim = IG::Entity::dimension;
 
@@ -280,6 +282,8 @@ public:
   void alpha_boundary (const IG& ig, const LFSU& lfsu_i, const X& x_i,
                                      const LFSV& lfsv_i, R& r_i) const
   {
+    assert_pointers();
+
     // Get entity dim
     const int dim = IG::Entity::dimension;
 
@@ -388,13 +392,16 @@ public:
       // Calculate diffusion coefficient from inside
       GradientU diff_coeff_i;
       hydr_disp_i.mv(normal_f,diff_coeff_i);
+      diff_coeff_i *= water_content_i;
+
+      // Use inside diffusion as face diffusion
+      const auto diff_coeff_f = diff_coeff_i*normal_f;
 
       // Inside unknown value
       RangeU u = x_i(lfsu_i,0);
 
       // Solute flux in normal direction w.r.t the intersection
-      auto soulte_flux_n = (g*water_flux_n)
-                          -(diff_coeff_i*water_content_i*(g-u)/distance);
+      auto soulte_flux_n = (g*water_flux_n)-(diff_coeff_f*(g-u)/distance);
 
       // Contribution to residual from Dirichlet boundary
       r_i.accumulate(lfsv_i, 0, soulte_flux_n*volume_f);
@@ -436,6 +443,7 @@ public:
    */
   void set_water_content(std::shared_ptr<GFWaterContent> gf_water_content)
   {
+      assert(gf_water_content);
       _gf_water_content = gf_water_content;
   }
 
@@ -446,11 +454,24 @@ public:
    */
   void set_water_flux(std::shared_ptr<GFWaterFlux> gf_water_flux)
   {
+      assert(gf_water_flux);
       _gf_water_flux = gf_water_flux;
   }
 
 private:
 
+  /// Assert that all shared pointers are assigned
+  /** \warning This uses the C `assert` macro and does nothing if compiled with
+   *           the `NDEBUG` definition (Release build).
+   */
+  void assert_pointers () const
+  {
+    assert(_param);
+    assert(_boundary);
+    assert(_gf_water_flux);
+    assert(_gf_water_content);
+  }
+
   const std::shared_ptr<const Parameter> _param;
   const std::shared_ptr<const Boundary> _boundary;
   std::shared_ptr<GFWaterFlux> _gf_water_flux;
@@ -528,6 +549,8 @@ public:
   void alpha_volume (const EG& eg, const LFSU& lfsu, const X& x, 
                                    const LFSV& lfsv, R& r) const
   {
+    assert(_gf_water_content);
+
     // Get local basis traits from local function space
     using LocalBasisTraitsU = typename LFSU::Traits::FiniteElementType::
                                             Traits::LocalBasisType::Traits;
@@ -572,6 +595,7 @@ public:
    */
   void set_water_content(std::shared_ptr<GFWaterContent> gf_water_content)
   {
+      assert(gf_water_content);
       _gf_water_content = gf_water_content;
   }