diff --git a/dumux/boxmodels/common/boxlocaljacobian.hh b/dumux/boxmodels/common/boxlocaljacobian.hh
index 41d4bfbe071fa4e2a8d4720d53ba5acd700c930d..9ccf1ffed81e285d28b74674d5e9d9fcff0eeab8 100644
--- a/dumux/boxmodels/common/boxlocaljacobian.hh
+++ b/dumux/boxmodels/common/boxlocaljacobian.hh
@@ -284,8 +284,18 @@ public:
     Scalar numericEpsilon(int scvIdx,
                           int pvIdx) const
     {
+        // define the base epsilon as the geometric mean of 1 and the
+        // resolution of the scalar type. E.g. for standard 64 bit
+        // floating point values, the resolution is about 10^-16 and
+        // the base epsilon is thus approximately 10^-8.
+        static const Scalar baseEps 
+            = Dumux::geometricMean<Scalar>(std::numeric_limits<Scalar>::epsilon(),
+                                           1.0);
+
+        // the epsilon value used for the numeric differentiation is
+        // now scaled by the absolute value of the primary variable...
         Scalar pv = this->curVolVars_[scvIdx].primaryVar(pvIdx);
-        return 1e-9*(std::abs(pv) + 1);
+        return baseEps*(std::abs(pv) + 1);
     }
 
 protected: