diff --git a/dumux/boxmodels/common/boxdarcyfluxvariables.hh b/dumux/boxmodels/common/boxdarcyfluxvariables.hh
index 953d7a3b62db147c1499a5cbfd613ff4ef0848a8..c7707293676d7f96660be763dcdc6bf063424d78 100644
--- a/dumux/boxmodels/common/boxdarcyfluxvariables.hh
+++ b/dumux/boxmodels/common/boxdarcyfluxvariables.hh
@@ -1,3 +1,78 @@
-#warning This file is deprecated. Include dumux/implicit/box/boxdarcyfluxvariables.hh instead.
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ *   See the file COPYING for full copying permissions.                      *
+ *                                                                           *
+ *   This program is free software: you can redistribute it and/or modify    *
+ *   it under the terms of the GNU General Public License as published by    *
+ *   the Free Software Foundation, either version 2 of the License, or       *
+ *   (at your option) any later version.                                     *
+ *                                                                           *
+ *   This program is distributed in the hope that it will be useful,         *
+ *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
+ *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
+ *   GNU General Public License for more details.                            *
+ *                                                                           *
+ *   You should have received a copy of the GNU General Public License       *
+ *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
+ *****************************************************************************/
+/*!
+ * \file
+ *
+ * \brief This file contains the data which is required to calculate
+ *        all fluxes of fluid phases over a face of a finite volume.
+ *
+ * This means pressure and temperature gradients, phase densities at
+ * the integration point, etc.
+ */
+#ifndef DUMUX_BOX_DARCY_FLUX_VARIABLES_HH
+#define DUMUX_BOX_DARCY_FLUX_VARIABLES_HH
 
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
+
+namespace Dumux
+{
+    
+/*!
+ * \ingroup BoxModel
+ * \ingroup BoxFluxVariables
+ * \brief Evaluates the normal component of the Darcy velocity 
+ * on a (sub)control volume face.
+ */
+template <class TypeTag>
+class BoxDarcyFluxVariables : public ImplicitDarcyFluxVariables<TypeTag>
+{
+    typedef ImplicitDarcyFluxVariables<TypeTag> ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GridView::template Codim<0>::Entity Element;
+    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+
+public:
+    /*
+     * \brief The constructor
+     *
+     * \param problem The problem
+     * \param element The finite element
+     * \param fvGeometry The finite-volume geometry in the box scheme
+     * \param faceIdx The local index of the SCV (sub-control-volume) face
+     * \param elemVolVars The volume variables of the current element
+     * \param onBoundary A boolean variable to specify whether the flux variables
+     * are calculated for interior SCV faces or boundary faces, default=false
+     */
+    DUNE_DEPRECATED_MSG("Use ImplicitDarcyFluxVariables from "
+                         "dumux/implicit/common/implicitdarcyfluxvariables.hh.")
+    BoxDarcyFluxVariables(const Problem &problem,
+                 const Element &element,
+                 const FVElementGeometry &fvGeometry,
+                 const int faceIdx,
+                 const ElementVolumeVariables &elemVolVars,
+                 const bool onBoundary = false)
+    : ParentType(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
+    {}
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/implicit/1p/1ppropertydefaults.hh b/dumux/implicit/1p/1ppropertydefaults.hh
index 34b3c254d4a93f329396d400c6e032c45759325e..7bacc08a40ce7e52aee7e14664cc36038abc9aa6 100644
--- a/dumux/implicit/1p/1ppropertydefaults.hh
+++ b/dumux/implicit/1p/1ppropertydefaults.hh
@@ -32,7 +32,7 @@
 #include "1pmodel.hh"
 #include "1plocalresidual.hh"
 #include "1pvolumevariables.hh"
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
 #include "1pindices.hh"
 
 #include <dumux/material/fluidsystems/gasphase.hh>
@@ -64,7 +64,7 @@ SET_TYPE_PROP(OneP, Model, OnePBoxModel<TypeTag>);
 SET_TYPE_PROP(OneP, VolumeVariables, OnePVolumeVariables<TypeTag>);
 
 //! the FluxVariables property
-SET_TYPE_PROP(OneP, FluxVariables, BoxDarcyFluxVariables<TypeTag>);
+SET_TYPE_PROP(OneP, FluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
 
 //! The indices required by the isothermal single-phase model
 SET_TYPE_PROP(OneP, Indices, OnePIndices);
diff --git a/dumux/implicit/1p/1pvolumevariables.hh b/dumux/implicit/1p/1pvolumevariables.hh
index c2e4588845a51278c23c9b897b2ffc1a8eb73285..01cbeab00b7809f6c3dee6c1fcdd9856c2ea75d3 100644
--- a/dumux/implicit/1p/1pvolumevariables.hh
+++ b/dumux/implicit/1p/1pvolumevariables.hh
@@ -145,7 +145,7 @@ public:
     /*!
      * \brief Returns the mobility.
      * 
-     * This function enables the use of BoxDarcyFluxVariables 
+     * This function enables the use of ImplicitDarcyFluxVariables 
      * with the 1p box model, ALTHOUGH the term mobility is
      * usually not employed in the one phase context. 
      *
diff --git a/dumux/implicit/2p/2ppropertydefaults.hh b/dumux/implicit/2p/2ppropertydefaults.hh
index ff139c3bf71d02da3336c09fcfccd285091a14c3..5fea8679838004dcbc5f7ae478eb2404ffc9db0a 100644
--- a/dumux/implicit/2p/2ppropertydefaults.hh
+++ b/dumux/implicit/2p/2ppropertydefaults.hh
@@ -30,7 +30,7 @@
 
 #include "2pmodel.hh"
 #include "2pindices.hh"
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
 #include "2pvolumevariables.hh"
 #include "2pproperties.hh"
 
@@ -68,7 +68,7 @@ SET_TYPE_PROP(TwoP, Model, TwoPModel<TypeTag>);
 SET_TYPE_PROP(TwoP, VolumeVariables, TwoPVolumeVariables<TypeTag>);
 
 //! the FluxVariables property
-SET_TYPE_PROP(TwoP, FluxVariables, BoxDarcyFluxVariables<TypeTag>);
+SET_TYPE_PROP(TwoP, FluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
 
 //! the upwind weight for the mass conservation equations.
 SET_SCALAR_PROP(TwoP, ImplicitMassUpwindWeight, 1.0);
diff --git a/dumux/implicit/2p2c/2p2cpropertydefaults.hh b/dumux/implicit/2p2c/2p2cpropertydefaults.hh
index 8a41d489879a0511f30b4be9af6e4ece8e49ed2e..7ecbdb6c4391c4f278e8a79ec790036e10b8a395 100644
--- a/dumux/implicit/2p2c/2p2cpropertydefaults.hh
+++ b/dumux/implicit/2p2c/2p2cpropertydefaults.hh
@@ -35,7 +35,7 @@
 #include "2p2cproperties.hh"
 #include "2p2cnewtoncontroller.hh"
 
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
 #include <dumux/material/spatialparams/boxspatialparams.hh>
 
 namespace Dumux
@@ -122,7 +122,7 @@ SET_TYPE_PROP(TwoPTwoC, VolumeVariables, TwoPTwoCVolumeVariables<TypeTag>);
 SET_TYPE_PROP(TwoPTwoC, FluxVariables, TwoPTwoCFluxVariables<TypeTag>);
 
 //! define the base flux variables to realize Darcy flow
-SET_TYPE_PROP(TwoPTwoC, BaseFluxVariables, BoxDarcyFluxVariables<TypeTag>);
+SET_TYPE_PROP(TwoPTwoC, BaseFluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
 
 //! the upwind weight for the mass conservation equations.
 SET_SCALAR_PROP(TwoPTwoC, ImplicitMassUpwindWeight, 1.0);
diff --git a/dumux/implicit/2pdfm/2pdfmfluxvariables.hh b/dumux/implicit/2pdfm/2pdfmfluxvariables.hh
index 0dba379ddb530eb924c549a1678f3b84391d77d1..958a06f317c20e6ddb06e846bf87dd9f18f866c7 100644
--- a/dumux/implicit/2pdfm/2pdfmfluxvariables.hh
+++ b/dumux/implicit/2pdfm/2pdfmfluxvariables.hh
@@ -27,7 +27,7 @@
 #include <dumux/common/math.hh>
 #include <dumux/common/parameters.hh>
 #include <dumux/common/valgrind.hh>
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
 
 #include "2pdfmproperties.hh"
 
@@ -45,7 +45,7 @@ namespace Dumux
  * the intergration point, etc.
  */
 template <class TypeTag>
-class TwoPDFMFluxVariables : public BoxDarcyFluxVariables<TypeTag>
+class TwoPDFMFluxVariables : public ImplicitDarcyFluxVariables<TypeTag>
 {
     typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
     typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
@@ -86,7 +86,7 @@ public:
                  int faceIdx,
                  const ElementVolumeVariables &elemVolVars,
                  const bool onBoundary = false)
-        : BoxDarcyFluxVariables<TypeTag>(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary), 
+        : ImplicitDarcyFluxVariables<TypeTag>(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary), 
           fvGeometry_(fvGeometry), faceIdx_(faceIdx), onBoundary_(onBoundary)
     {
         faceSCV_ = &this->face();
diff --git a/dumux/implicit/2pdfm/2pdfmpropertydefaults.hh b/dumux/implicit/2pdfm/2pdfmpropertydefaults.hh
index 1af745335d0c1638e22aed0a25dc58096b2566c7..0ea4c8d9bdeaa7ecbbc493584bf1dea1f214ba62 100644
--- a/dumux/implicit/2pdfm/2pdfmpropertydefaults.hh
+++ b/dumux/implicit/2pdfm/2pdfmpropertydefaults.hh
@@ -27,7 +27,7 @@
 #ifndef DUMUX_BOXMODELS_2PDFM_PROPERTY_DEFAULTS_HH
 #define DUMUX_BOXMODELS_2PDFM_PROPERTY_DEFAULTS_HH
 
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
 #include <dumux/material/components/nullcomponent.hh>
 #include <dumux/material/fluidstates/immisciblefluidstate.hh>
 #include <dumux/material/fluidsystems/2pimmisciblefluidsystem.hh>
diff --git a/dumux/implicit/2pni/2pnifluxvariables.hh b/dumux/implicit/2pni/2pnifluxvariables.hh
index 2d331d6a89e7d3d7b7b53d4ad0fe70d90e4c6ef5..ca7496b340d5d17744e2901b3ee097361ddae110 100644
--- a/dumux/implicit/2pni/2pnifluxvariables.hh
+++ b/dumux/implicit/2pni/2pnifluxvariables.hh
@@ -29,7 +29,7 @@
 #define DUMUX_2PNI_FLUX_VARIABLES_HH
 
 #include <dumux/common/math.hh>
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
 
 namespace Dumux
 {
@@ -45,9 +45,9 @@ namespace Dumux
  * the integration point, etc.
  */
 template <class TypeTag>
-class TwoPNIFluxVariables : public BoxDarcyFluxVariables<TypeTag>
+class TwoPNIFluxVariables : public ImplicitDarcyFluxVariables<TypeTag>
 {
-    typedef BoxDarcyFluxVariables<TypeTag> ParentType;
+    typedef ImplicitDarcyFluxVariables<TypeTag> ParentType;
 
     typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
     typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
diff --git a/dumux/implicit/3p3c/3p3cpropertydefaults.hh b/dumux/implicit/3p3c/3p3cpropertydefaults.hh
index a644c778c93a1d9fef2dfe6df7b2b7fe899fb891..4309856760af22d5e5169cc483b773c184b6292d 100644
--- a/dumux/implicit/3p3c/3p3cpropertydefaults.hh
+++ b/dumux/implicit/3p3c/3p3cpropertydefaults.hh
@@ -38,7 +38,7 @@
 #include "3p3cnewtoncontroller.hh"
 // #include "3p3cboundaryvariables.hh"
 
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
 #include <dumux/material/spatialparams/boxspatialparams.hh>
 
 namespace Dumux
@@ -108,7 +108,7 @@ SET_TYPE_PROP(ThreePThreeC, VolumeVariables, ThreePThreeCVolumeVariables<TypeTag
 SET_TYPE_PROP(ThreePThreeC, FluxVariables, ThreePThreeCFluxVariables<TypeTag>);
 
 //! define the base flux variables to realize Darcy flow
-SET_TYPE_PROP(ThreePThreeC, BaseFluxVariables, BoxDarcyFluxVariables<TypeTag>);
+SET_TYPE_PROP(ThreePThreeC, BaseFluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
 
 //! the upwind factor for the mobility.
 SET_SCALAR_PROP(ThreePThreeC, ImplicitMassUpwindWeight, 1.0);
diff --git a/dumux/implicit/box/boxdarcyfluxvariables.hh b/dumux/implicit/box/boxdarcyfluxvariables.hh
deleted file mode 100644
index e950b5323ec5eea0ed052e5473f033e73b590843..0000000000000000000000000000000000000000
--- a/dumux/implicit/box/boxdarcyfluxvariables.hh
+++ /dev/null
@@ -1,324 +0,0 @@
-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- *   See the file COPYING for full copying permissions.                      *
- *                                                                           *
- *   This program is free software: you can redistribute it and/or modify    *
- *   it under the terms of the GNU General Public License as published by    *
- *   the Free Software Foundation, either version 2 of the License, or       *
- *   (at your option) any later version.                                     *
- *                                                                           *
- *   This program is distributed in the hope that it will be useful,         *
- *   but WITHOUT ANY WARRANTY; without even the implied warranty of          *
- *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the            *
- *   GNU General Public License for more details.                            *
- *                                                                           *
- *   You should have received a copy of the GNU General Public License       *
- *   along with this program.  If not, see <http://www.gnu.org/licenses/>.   *
- *****************************************************************************/
-/*!
- * \file
- *
- * \brief This file contains the data which is required to calculate
- *        all fluxes of fluid phases over a face of a finite volume.
- *
- * This means pressure and temperature gradients, phase densities at
- * the integration point, etc.
- */
-#ifndef DUMUX_BOX_DARCY_FLUX_VARIABLES_HH
-#define DUMUX_BOX_DARCY_FLUX_VARIABLES_HH
-
-#include "boxproperties.hh"
-
-#include <dumux/common/parameters.hh>
-#include <dumux/common/math.hh>
-
-namespace Dumux
-{
-    
-namespace Properties
-{
-// forward declaration of properties 
-NEW_PROP_TAG(ImplicitMobilityUpwindWeight);
-NEW_PROP_TAG(SpatialParams);
-NEW_PROP_TAG(NumPhases);
-NEW_PROP_TAG(ProblemEnableGravity);
-}   
-
-/*!
- * \ingroup BoxModel
- * \ingroup BoxFluxVariables
- * \brief Evaluates the normal component of the Darcy velocity 
- * on a (sub)control volume face.
- */
-template <class TypeTag>
-class BoxDarcyFluxVariables
-{
-    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-    typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
-    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
-    
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-    typedef typename GridView::template Codim<0>::Entity Element;
-
-    enum { dim = GridView::dimension} ;
-    enum { dimWorld = GridView::dimensionworld} ;
-    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef Dune::FieldMatrix<Scalar, dimWorld, dimWorld> Tensor;
-    typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
-
-    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-    typedef typename FVElementGeometry::SubControlVolumeFace SCVFace;
-
-public:
-    /*
-     * \brief The constructor
-     *
-     * \param problem The problem
-     * \param element The finite element
-     * \param fvGeometry The finite-volume geometry in the box scheme
-     * \param faceIdx The local index of the SCV (sub-control-volume) face
-     * \param elemVolVars The volume variables of the current element
-     * \param onBoundary A boolean variable to specify whether the flux variables
-     * are calculated for interior SCV faces or boundary faces, default=false
-     */
-    BoxDarcyFluxVariables(const Problem &problem,
-                 const Element &element,
-                 const FVElementGeometry &fvGeometry,
-                 const int faceIdx,
-                 const ElementVolumeVariables &elemVolVars,
-                 const bool onBoundary = false)
-    : fvGeometry_(fvGeometry), faceIdx_(faceIdx), onBoundary_(onBoundary)
-    {
-        mobilityUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MobilityUpwindWeight);
-        calculateGradients_(problem, element, elemVolVars);
-        calculateNormalVelocity_(problem, element, elemVolVars);
-    }
-
-public:
-    /*!
-     * \brief Return the volumetric flux over a face of a given phase.
-     *
-     *        This is the calculated velocity multiplied by the unit normal
-     *        and the area of the face.
-     *        face().normal
-     *        has already the magnitude of the area.
-     *
-     * \param phaseIdx index of the phase
-     */
-    Scalar volumeFlux(const unsigned int phaseIdx) const
-    { return volumeFlux_[phaseIdx]; }
-    
-    /*!
-     * \brief Return the velocity of a given phase.
-     *
-     *        This is the full velocity vector on the
-     *        face (without being multiplied with normal).
-     *
-     * \param phaseIdx index of the phase
-     */
-    DimVector velocity(const unsigned int phaseIdx) const
-    { return velocity_[phaseIdx] ; }
-
-    /*!
-     * \brief Return intrinsic permeability multiplied with potential
-     *        gradient multiplied with normal.
-     *        I.e. everything that does not need upwind decisions.
-     *
-     * \param phaseIdx index of the phase
-     */
-    Scalar kGradPNormal(const unsigned int phaseIdx) const
-    { return kGradPNormal_[phaseIdx] ; }
-
-    /*!
-     * \brief Return the local index of the downstream control volume
-     *        for a given phase.
-     *
-     * \param phaseIdx index of the phase
-     */
-    const unsigned int downstreamIdx(const unsigned phaseIdx) const
-    { return downstreamIdx_[phaseIdx]; }
-    
-    /*!
-     * \brief Return the local index of the upstream control volume
-     *        for a given phase.
-     *
-     * \param phaseIdx index of the phase
-     */
-    const unsigned int upstreamIdx(const unsigned phaseIdx) const
-    { return upstreamIdx_[phaseIdx]; }
-
-    /*!
-     * \brief Return the SCV (sub-control-volume) face. This may be either
-     *        a face within the element or a face on the element boundary,
-     *        depending on the value of onBoundary_.
-     */
-    const SCVFace &face() const
-    {
-        if (onBoundary_)
-            return fvGeometry_.boundaryFace[faceIdx_];
-        else
-            return fvGeometry_.subContVolFace[faceIdx_];
-    }
-
-protected:
-    
-    /*
-     * \brief Calculation of the potential gradients
-     *
-     * \param problem The problem
-     * \param element The finite element
-     * \param elemVolVars The volume variables of the current element
-     * are calculated for interior SCV faces or boundary faces, default=false
-     */
-    void calculateGradients_(const Problem &problem,
-                             const Element &element,
-                             const ElementVolumeVariables &elemVolVars)
-    {
-        // loop over all phases
-        for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
-        {
-            gradPotential_[phaseIdx]= 0.0 ;
-            for (int idx = 0;
-                 idx < fvGeometry_.numFAP;
-                 idx++) // loop over adjacent vertices
-            {
-                // FE gradient at vertex idx
-                const DimVector &feGrad = face().grad[idx];
-
-                // index for the element volume variables
-                int volVarsIdx = face().fapIndices[idx];
-
-                // the pressure gradient
-                DimVector tmp(feGrad);
-                tmp *= elemVolVars[volVarsIdx].fluidState().pressure(phaseIdx);
-                gradPotential_[phaseIdx] += tmp;
-            }
-
-            // correct the pressure gradient by the gravitational acceleration
-            if (GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity))
-            {
-                // estimate the gravitational acceleration at a given SCV face
-                // using the arithmetic mean
-                DimVector g(problem.boxGravity(element, fvGeometry_, face().i));
-                g += problem.boxGravity(element, fvGeometry_, face().j);
-                g /= 2;
-
-                // calculate the phase density at the integration point. we
-                // only do this if the wetting phase is present in both cells
-                Scalar SI = elemVolVars[face().i].fluidState().saturation(phaseIdx);
-                Scalar SJ = elemVolVars[face().j].fluidState().saturation(phaseIdx);
-                Scalar rhoI = elemVolVars[face().i].fluidState().density(phaseIdx);
-                Scalar rhoJ = elemVolVars[face().j].fluidState().density(phaseIdx);
-                Scalar fI = std::max(0.0, std::min(SI/1e-5, 0.5));
-                Scalar fJ = std::max(0.0, std::min(SJ/1e-5, 0.5));
-                if (fI + fJ == 0)
-                    // doesn't matter because no wetting phase is present in
-                    // both cells!
-                    fI = fJ = 0.5;
-                const Scalar density = (fI*rhoI + fJ*rhoJ)/(fI + fJ);
-
-                // make gravity acceleration a force
-                DimVector f(g);
-                f *= density;
-
-                // calculate the final potential gradient
-                gradPotential_[phaseIdx] -= f;
-            } // gravity
-        } // loop over all phases
-     }
-
-    /*
-     * \brief Actual calculation of the normal Darcy velocities.
-     *
-     * \param problem The problem
-     * \param element The finite element
-     * \param elemVolVars The volume variables of the current element
-     */
-    void calculateNormalVelocity_(const Problem &problem,
-                                  const Element &element,
-                                  const ElementVolumeVariables &elemVolVars)
-    {
-        // calculate the mean intrinsic permeability
-        const SpatialParams &spatialParams = problem.spatialParams();
-        Tensor K;
-        if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
-        {
-            spatialParams.meanK(K,
-                                spatialParams.intrinsicPermeability(element,
-                                                                    fvGeometry_,
-                                                                    face().i),
-                                spatialParams.intrinsicPermeability(element,
-                                                                    fvGeometry_,
-                                                                    face().j));
-        }
-        else
-        {
-            spatialParams.meanK(K,
-                                spatialParams.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().i]),
-                                spatialParams.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().j]));
-        }
-        
-        // loop over all phases
-        for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
-        {
-            // calculate the flux in the normal direction of the
-            // current sub control volume face:
-            //
-            // v = - (K_f grad phi) * n
-            // with K_f = rho g / mu K
-            //
-            // Mind, that the normal has the length of it's area.
-            // This means that we are actually calculating
-            //  Q = - (K grad phi) dot n /|n| * A
-
-
-            K.mv(gradPotential_[phaseIdx], kGradP_[phaseIdx]);
-            kGradPNormal_[phaseIdx] = kGradP_[phaseIdx]*face().normal;
-
-            // determine the upwind direction
-            if (kGradPNormal_[phaseIdx] < 0)
-            {
-                upstreamIdx_[phaseIdx] = face().i;
-                downstreamIdx_[phaseIdx] = face().j;
-            }
-            else
-            {
-                upstreamIdx_[phaseIdx] = face().j;
-                downstreamIdx_[phaseIdx] = face().i;
-            }
-            
-            // obtain the upwind volume variables
-            const VolumeVariables& upVolVars = elemVolVars[ upstreamIdx(phaseIdx) ];
-            const VolumeVariables& downVolVars = elemVolVars[ downstreamIdx(phaseIdx) ];
-
-            // the minus comes from the Darcy relation which states that
-            // the flux is from high to low potentials.
-            // set the velocity
-            velocity_[phaseIdx] = kGradP_[phaseIdx];
-            velocity_[phaseIdx] *= - ( mobilityUpwindWeight_*upVolVars.mobility(phaseIdx)
-                    + (1.0 - mobilityUpwindWeight_)*downVolVars.mobility(phaseIdx)) ;
-
-            // set the volume flux
-            volumeFlux_[phaseIdx] = velocity_[phaseIdx] * face().normal ;
-        }// loop all phases
-    }
-
-    const FVElementGeometry &fvGeometry_;   //!< Information about the geometry of discretization
-    const unsigned int faceIdx_;            //!< The index of the sub control volume face
-    const bool      onBoundary_;                //!< Specifying whether we are currently on the boundary of the simulation domain
-    unsigned int    upstreamIdx_[numPhases] , downstreamIdx_[numPhases]; //!< local index of the upstream / downstream vertex
-    Scalar          volumeFlux_[numPhases] ;    //!< Velocity multiplied with normal (magnitude=area)
-    DimVector       velocity_[numPhases] ;      //!< The velocity as determined by Darcy's law or by the Forchheimer relation
-    Scalar          kGradPNormal_[numPhases] ;  //!< Permeability multiplied with gradient in potential, multiplied with normal (magnitude=area)
-    DimVector       kGradP_[numPhases] ; //!< Permeability multiplied with gradient in potential
-    DimVector       gradPotential_[numPhases] ; //!< Gradient of potential, which drives flow
-    Scalar          mobilityUpwindWeight_;      //!< Upwind weight for mobility. Set to one for full upstream weighting
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/implicit/box/boxforchheimerfluxvariables.hh b/dumux/implicit/box/boxforchheimerfluxvariables.hh
index 78d008f20316903f1c6235e14cbc7bea4d4c3fd5..98b4334197c431b25e449c296d20cf84fcb90216 100644
--- a/dumux/implicit/box/boxforchheimerfluxvariables.hh
+++ b/dumux/implicit/box/boxforchheimerfluxvariables.hh
@@ -31,7 +31,7 @@
 
 #include <dumux/common/parameters.hh>
 #include <dumux/common/math.hh>
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
 
 namespace Dumux
 {
@@ -76,7 +76,7 @@ NEW_PROP_TAG(ProblemEnableGravity);
  */
 template <class TypeTag>
 class BoxForchheimerFluxVariables
-    : public BoxDarcyFluxVariables<TypeTag>
+    : public ImplicitDarcyFluxVariables<TypeTag>
 {
     typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
     typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
@@ -113,7 +113,7 @@ public:
                  const unsigned int faceIdx,
                  const ElementVolumeVariables &elemVolVars,
                  const bool onBoundary = false)
-        :   BoxDarcyFluxVariables<TypeTag>(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
+        :   ImplicitDarcyFluxVariables<TypeTag>(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
     {
         calculateNormalVelocity_(problem, element, elemVolVars);
     }
diff --git a/dumux/implicit/common/implicitdarcyfluxvariables.hh b/dumux/implicit/common/implicitdarcyfluxvariables.hh
index 6e5d3df449f65a22a12966497c6a5a2c0a6e800c..d529fce5e4bc94b2b4115e52efc3f38666543bfe 100644
--- a/dumux/implicit/common/implicitdarcyfluxvariables.hh
+++ b/dumux/implicit/common/implicitdarcyfluxvariables.hh
@@ -25,10 +25,10 @@
  * This means pressure and temperature gradients, phase densities at
  * the integration point, etc.
  */
-#ifndef DUMUX_BOX_DARCY_FLUX_VARIABLES_HH
-#define DUMUX_BOX_DARCY_FLUX_VARIABLES_HH
+#ifndef DUMUX_IMPLICIT_DARCY_FLUX_VARIABLES_HH
+#define DUMUX_IMPLICIT_DARCY_FLUX_VARIABLES_HH
 
-#include "boxproperties.hh"
+#include "implicitproperties.hh"
 
 #include <dumux/common/parameters.hh>
 #include <dumux/common/math.hh>
@@ -46,13 +46,13 @@ NEW_PROP_TAG(ProblemEnableGravity);
 }   
 
 /*!
- * \ingroup BoxModel
- * \ingroup BoxFluxVariables
+ * \ingroup ImplicitModel
+ * \ingroup ImplicitFluxVariables
  * \brief Evaluates the normal component of the Darcy velocity 
  * on a (sub)control volume face.
  */
 template <class TypeTag>
-class BoxDarcyFluxVariables
+class ImplicitDarcyFluxVariables
 {
     typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
     typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
@@ -85,13 +85,13 @@ public:
      * \param onBoundary A boolean variable to specify whether the flux variables
      * are calculated for interior SCV faces or boundary faces, default=false
      */
-    BoxDarcyFluxVariables(const Problem &problem,
+    ImplicitDarcyFluxVariables(const Problem &problem,
                  const Element &element,
                  const FVElementGeometry &fvGeometry,
                  const int faceIdx,
                  const ElementVolumeVariables &elemVolVars,
                  const bool onBoundary = false)
-        : fvGeometry_(fvGeometry), faceIdx_(faceIdx), onBoundary_(onBoundary)
+    : fvGeometry_(fvGeometry), faceIdx_(faceIdx), onBoundary_(onBoundary)
     {
         mobilityUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MobilityUpwindWeight);
         calculateGradients_(problem, element, elemVolVars);
@@ -165,17 +165,7 @@ public:
     }
 
 protected:
-    const FVElementGeometry &fvGeometry_;   //!< Information about the geometry of discretization
-    const unsigned int faceIdx_;            //!< The index of the sub control volume face
-    const bool      onBoundary_;                //!< Specifying whether we are currently on the boundary of the simulation domain
-    unsigned int    upstreamIdx_[numPhases] , downstreamIdx_[numPhases]; //!< local index of the upstream / downstream vertex
-    Scalar          volumeFlux_[numPhases] ;    //!< Velocity multiplied with normal (magnitude=area)
-    DimVector       velocity_[numPhases] ;      //!< The velocity as determined by Darcy's law or by the Forchheimer relation
-    Scalar          kGradPNormal_[numPhases] ;  //!< Permeability multiplied with gradient in potential, multiplied with normal (magnitude=area)
-    DimVector       kGradP_[numPhases] ; //!< Permeability multiplied with gradient in potential
-    DimVector       gradPotential_[numPhases] ; //!< Gradient of potential, which drives flow
-    Scalar          mobilityUpwindWeight_;      //!< Upwind weight for mobility. Set to one for full upstream weighting
-
+    
     /*
      * \brief Calculation of the potential gradients
      *
@@ -255,13 +245,23 @@ protected:
         // calculate the mean intrinsic permeability
         const SpatialParams &spatialParams = problem.spatialParams();
         Tensor K;
-        spatialParams.meanK(K,
-                            spatialParams.intrinsicPermeability(element,
-                                                                fvGeometry_,
-                                                                face().i),
-                            spatialParams.intrinsicPermeability(element,
-                                                                fvGeometry_,
-                                                                face().j));
+        if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
+        {
+            spatialParams.meanK(K,
+                                spatialParams.intrinsicPermeability(element,
+                                                                    fvGeometry_,
+                                                                    face().i),
+                                spatialParams.intrinsicPermeability(element,
+                                                                    fvGeometry_,
+                                                                    face().j));
+        }
+        else
+        {
+            spatialParams.meanK(K,
+                                spatialParams.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().i]),
+                                spatialParams.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().j]));
+        }
+        
         // loop over all phases
         for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
         {
@@ -306,6 +306,17 @@ protected:
             volumeFlux_[phaseIdx] = velocity_[phaseIdx] * face().normal ;
         }// loop all phases
     }
+
+    const FVElementGeometry &fvGeometry_;   //!< Information about the geometry of discretization
+    const unsigned int faceIdx_;            //!< The index of the sub control volume face
+    const bool      onBoundary_;                //!< Specifying whether we are currently on the boundary of the simulation domain
+    unsigned int    upstreamIdx_[numPhases] , downstreamIdx_[numPhases]; //!< local index of the upstream / downstream vertex
+    Scalar          volumeFlux_[numPhases] ;    //!< Velocity multiplied with normal (magnitude=area)
+    DimVector       velocity_[numPhases] ;      //!< The velocity as determined by Darcy's law or by the Forchheimer relation
+    Scalar          kGradPNormal_[numPhases] ;  //!< Permeability multiplied with gradient in potential, multiplied with normal (magnitude=area)
+    DimVector       kGradP_[numPhases] ; //!< Permeability multiplied with gradient in potential
+    DimVector       gradPotential_[numPhases] ; //!< Gradient of potential, which drives flow
+    Scalar          mobilityUpwindWeight_;      //!< Upwind weight for mobility. Set to one for full upstream weighting
 };
 
 } // end namespace
diff --git a/dumux/implicit/common/implicitforchheimerfluxvariables.hh b/dumux/implicit/common/implicitforchheimerfluxvariables.hh
index ec3117cb08f9215b4d3048aaec4f665abbb11de9..98b4334197c431b25e449c296d20cf84fcb90216 100644
--- a/dumux/implicit/common/implicitforchheimerfluxvariables.hh
+++ b/dumux/implicit/common/implicitforchheimerfluxvariables.hh
@@ -31,7 +31,7 @@
 
 #include <dumux/common/parameters.hh>
 #include <dumux/common/math.hh>
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
 
 namespace Dumux
 {
@@ -76,7 +76,7 @@ NEW_PROP_TAG(ProblemEnableGravity);
  */
 template <class TypeTag>
 class BoxForchheimerFluxVariables
-    : public BoxDarcyFluxVariables<TypeTag>
+    : public ImplicitDarcyFluxVariables<TypeTag>
 {
     typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
     typedef typename GET_PROP_TYPE(TypeTag, SpatialParams) SpatialParams;
@@ -113,7 +113,7 @@ public:
                  const unsigned int faceIdx,
                  const ElementVolumeVariables &elemVolVars,
                  const bool onBoundary = false)
-        :   BoxDarcyFluxVariables<TypeTag>(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
+        :   ImplicitDarcyFluxVariables<TypeTag>(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
     {
         calculateNormalVelocity_(problem, element, elemVolVars);
     }
@@ -133,14 +133,23 @@ protected:
         // calculate the mean intrinsic permeability
         const SpatialParams &spatialParams = problem.spatialParams();
         Tensor K;
-        spatialParams.meanK(K,
-                            spatialParams.intrinsicPermeability(element,
-                                                                this->fvGeometry_,
-                                                                this->face().i),
-                            spatialParams.intrinsicPermeability(element,
-                                                                this->fvGeometry_,
-                                                                this->face().j));
-
+        if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
+        {
+            spatialParams.meanK(K,
+                                spatialParams.intrinsicPermeability(element,
+                                                                    fvGeometry_,
+                                                                    face().i),
+                                spatialParams.intrinsicPermeability(element,
+                                                                    fvGeometry_,
+                                                                    face().j));
+        }
+        else
+        {
+            spatialParams.meanK(K,
+                                spatialParams.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().i]),
+                                spatialParams.elementIntrinsicPermeability(*fvGeometry_.neighbors[face().j]));
+        }
+        
         // obtain the Forchheimer coefficient from the spatial parameters
         const Scalar forchCoeff = spatialParams.forchCoeff(element,
                                                           this->fvGeometry_,
diff --git a/dumux/implicit/mpnc/mpncfluxvariables.hh b/dumux/implicit/mpnc/mpncfluxvariables.hh
index b557dcd1437c733b7282f8aca85afeea21355685..cfc9a992b2bc5ab187e71d6202713605e55b8df7 100644
--- a/dumux/implicit/mpnc/mpncfluxvariables.hh
+++ b/dumux/implicit/mpnc/mpncfluxvariables.hh
@@ -32,7 +32,7 @@
 
 #include "diffusion/fluxvariables.hh"
 #include "energy/mpncfluxvariablesenergy.hh"
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
 #include <dumux/implicit/box/boxforchheimerfluxvariables.hh>
 
 namespace Dumux
diff --git a/dumux/implicit/mpnc/mpncmodel.hh b/dumux/implicit/mpnc/mpncmodel.hh
index 18f9cd2e71d43790a9f02a326b3dc7ef2eeb42a6..83ac43108407465facef8e2be6e519c42650fcea 100644
--- a/dumux/implicit/mpnc/mpncmodel.hh
+++ b/dumux/implicit/mpnc/mpncmodel.hh
@@ -38,7 +38,7 @@ namespace Dumux
  * denoted by the upper index \f$\kappa \in \{ 1, \dots, N \} \f$.
  *
  * The momentum approximation can be selected via "BaseFluxVariables":
- * Darcy (BoxDarcyFluxVariables) and Forchheimer (BoxForchheimerFluxVariables)
+ * Darcy (ImplicitDarcyFluxVariables) and Forchheimer (BoxForchheimerFluxVariables)
  * relations are available for all Box models.
  *
  * By inserting this into the equations for the conservation of the
diff --git a/dumux/implicit/mpnc/mpncpropertydefaults.hh b/dumux/implicit/mpnc/mpncpropertydefaults.hh
index 04c5948613b6a9243f97f1c7dcbcee0b28297c49..9a910ea32da1ef973a6346a289615b657a870da0 100644
--- a/dumux/implicit/mpnc/mpncpropertydefaults.hh
+++ b/dumux/implicit/mpnc/mpncpropertydefaults.hh
@@ -173,7 +173,7 @@ SET_TYPE_PROP(MPNC, VolumeVariables, MPNCVolumeVariables<TypeTag>);
 SET_TYPE_PROP(MPNC, FluxVariables, MPNCFluxVariables<TypeTag>);
 
 //! the Base of the Fluxvariables property (Darcy / Forchheimer)
-SET_TYPE_PROP(MPNC, BaseFluxVariables, BoxDarcyFluxVariables<TypeTag>);
+SET_TYPE_PROP(MPNC, BaseFluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
 
 //! truncate the newton update for the first few Newton iterations of a time step
 SET_BOOL_PROP(MPNC, NewtonEnableChop, true);
diff --git a/dumux/implicit/richards/richardspropertydefaults.hh b/dumux/implicit/richards/richardspropertydefaults.hh
index 2e312d08051bbb19dec0631e18087b387804a5af..2766eb2ffff0247b15ffb960c73fd62364772930 100644
--- a/dumux/implicit/richards/richardspropertydefaults.hh
+++ b/dumux/implicit/richards/richardspropertydefaults.hh
@@ -31,7 +31,7 @@
 #include "richardsmodel.hh"
 #include "richardsproblem.hh"
 #include "richardsindices.hh"
-#include <dumux/implicit/box/boxdarcyfluxvariables.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
 #include "richardsvolumevariables.hh"
 #include "richardsproperties.hh"
 #include "richardsnewtoncontroller.hh"
@@ -65,7 +65,7 @@ SET_TYPE_PROP(Richards, Model, RichardsModel<TypeTag>);
 SET_TYPE_PROP(Richards, VolumeVariables, RichardsVolumeVariables<TypeTag>);
 
 //! The class for the quantities required for the flux calculation
-SET_TYPE_PROP(Richards, FluxVariables, BoxDarcyFluxVariables<TypeTag>);
+SET_TYPE_PROP(Richards, FluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
 
 //! The class of the newton controller
 SET_TYPE_PROP(Richards, NewtonController, RichardsNewtonController<TypeTag>);