diff --git a/dumux/implicit/2p/2pproperties.hh b/dumux/implicit/2p/2pproperties.hh
index f3e78f16b84e7651356e48651f69444d7e10c5e0..d650e9607ecface7bc08ae06c85d4a671c178619 100644
--- a/dumux/implicit/2p/2pproperties.hh
+++ b/dumux/implicit/2p/2pproperties.hh
@@ -54,13 +54,10 @@ NEW_TYPE_TAG(TwoP);
 NEW_TYPE_TAG(BoxTwoP, INHERITS_FROM(BoxModel, TwoP));
 NEW_TYPE_TAG(CCTwoP, INHERITS_FROM(CCModel, TwoP));
 
-// TODO: remove the following ifndef after the old 2pni model has been removed
-#ifndef DUMUX_2PNI_PROPERTIES_HH
 //! The type tags for the corresponding non-isothermal problems
 NEW_TYPE_TAG(TwoPNI, INHERITS_FROM(TwoP, NonIsothermal));
 NEW_TYPE_TAG(BoxTwoPNI, INHERITS_FROM(BoxModel, TwoPNI));
 NEW_TYPE_TAG(CCTwoPNI, INHERITS_FROM(CCModel, TwoPNI));
-#endif
 
 //////////////////////////////////////////////////////////////////
 // Property tags
diff --git a/dumux/implicit/2p2c/2p2cproperties.hh b/dumux/implicit/2p2c/2p2cproperties.hh
index 69e8f46c129f09cf563f52425ae04ce95b1185d3..da4f9aca3ff31a63023904e890a9b1a8da18c8ca 100644
--- a/dumux/implicit/2p2c/2p2cproperties.hh
+++ b/dumux/implicit/2p2c/2p2cproperties.hh
@@ -47,13 +47,10 @@ NEW_TYPE_TAG(TwoPTwoC);
 NEW_TYPE_TAG(BoxTwoPTwoC, INHERITS_FROM(BoxModel, TwoPTwoC));
 NEW_TYPE_TAG(CCTwoPTwoC, INHERITS_FROM(CCModel, TwoPTwoC));
 
-// TODO: remove the following ifndef after the old 2p2cni model has been removed
-#ifndef DUMUX_2P2CNI_PROPERTIES_HH
 //! The type tags for the corresponding non-isothermal problems
 NEW_TYPE_TAG(TwoPTwoCNI, INHERITS_FROM(TwoPTwoC, NonIsothermal));
 NEW_TYPE_TAG(BoxTwoPTwoCNI, INHERITS_FROM(BoxModel, TwoPTwoCNI));
 NEW_TYPE_TAG(CCTwoPTwoCNI, INHERITS_FROM(CCModel, TwoPTwoCNI));
-#endif
 
 //////////////////////////////////////////////////////////////////
 // Property tags
diff --git a/dumux/implicit/2p2cni/2p2cnifluxvariables.hh b/dumux/implicit/2p2cni/2p2cnifluxvariables.hh
deleted file mode 100644
index 8209c56e9e07efb65722b6974d49481f88cbc333..0000000000000000000000000000000000000000
--- a/dumux/implicit/2p2cni/2p2cnifluxvariables.hh
+++ /dev/null
@@ -1,205 +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 data which is required to calculate
- *        the heat fluxes over a face of a finite volume.
- *
- * This means temperature gradients and the normal matrix
- * heat flux.
- */
-#ifndef DUMUX_2P2CNI_FLUX_VARIABLES_HH
-#define DUMUX_2P2CNI_FLUX_VARIABLES_HH
-
-#include "2p2cniproperties.hh"
-#include <dumux/common/math.hh>
-#include <dumux/implicit/2p2c/2p2cfluxvariables.hh>
-
-namespace Dumux
-{
-
-/*!
- * \ingroup TwoPTwoCNIModel
- * \ingroup ImplicitFluxVariables
- * \brief This template class contains data which is required to
- *        calculate the heat fluxes over a face of a finite
- *        volume for the non-isothermal two-phase two-component model.
- *        The mass fluxes are computed in the parent class.
- *
- * This means temperature gradients and the normal matrix
- * heat flux.
- */
-template <class TypeTag>
-class TwoPTwoCNIFluxVariables : public TwoPTwoCFluxVariables<TypeTag>
-{
-    typedef TwoPTwoCFluxVariables<TypeTag> ParentType;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-    typedef typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel) ThermalConductivityModel;
-
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-    typedef typename GridView::template Codim<0>::Entity Element;
-    enum { dimWorld = GridView::dimensionworld };
-    enum { dim = GridView::dimension };
-    typedef Dune::FieldVector<Scalar, dim> DimVector;
-    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-    enum {
-        wPhaseIdx = Indices::wPhaseIdx,
-        nPhaseIdx = Indices::nPhaseIdx
-    };
-
-public:
-    /*!
-     * \brief The constructor
-     *
-     * \param problem The problem
-     * \param element The finite element
-     * \param fvGeometry The finite-volume geometry in the fully implicit scheme
-     * \param faceIdx The local index of the sub-control-volume face
-     * \param elemVolVars The volume variables of the current element
-     * \param onBoundary Distinguishes if we are on a sub-control-volume face or on a boundary face
-     */
-    TwoPTwoCNIFluxVariables(const Problem &problem,
-                            const Element &element,
-                            const FVElementGeometry &fvGeometry,
-                            const int faceIdx,
-                            const ElementVolumeVariables &elemVolVars,
-                            bool onBoundary = false)
-        : ParentType(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
-    {
-        faceIdx_ = faceIdx;
-
-        calculateValues_(problem, element, elemVolVars);
-    }
-
-    /*!
-     * \brief The total heat flux \f$\mathrm{[J/s]}\f$ due to heat conduction
-     *        of the rock matrix over the sub-control volume face in
-     *        direction of the face normal.
-     */
-    Scalar normalMatrixHeatFlux() const
-    { return normalMatrixHeatFlux_; }
-
-    /*!
-     * \brief The local temperature gradient at the IP of the considered scv face.
-     */
-    GlobalPosition temperatureGradient() const
-    { return temperatureGrad_; }
-
-    /*!
-     * \brief The harmonically averaged effective thermal conductivity.
-     */
-    Scalar effThermalConductivity() const
-    { return lambdaEff_; }
-
-protected:
-    void calculateValues_(const Problem &problem,
-                          const Element &element,
-                          const ElementVolumeVariables &elemVolVars)
-    {
-        // calculate temperature gradient using finite element
-        // gradients
-        temperatureGrad_ = 0;
-        GlobalPosition tmp(0.0);
-        for (unsigned int idx = 0; idx < this->face().numFap; idx++)
-        {
-            tmp = this->face().grad[idx];
-
-            // index for the element volume variables
-            int volVarsIdx = this->face().fapIndices[idx];
-
-            tmp *= elemVolVars[volVarsIdx].temperature();
-            temperatureGrad_ += tmp;
-        }
-
-        lambdaEff_ = 0;
-        calculateEffThermalConductivity_(problem, element, elemVolVars);
-
-        // project the heat flux vector on the face's normal vector
-        normalMatrixHeatFlux_ = temperatureGrad_* this->face().normal;
-        normalMatrixHeatFlux_ *= -lambdaEff_;
-    }
-
-    void calculateEffThermalConductivity_(const Problem &problem,
-                                          const Element &element,
-                                          const ElementVolumeVariables &elemVolVars)
-    {
-        const unsigned i = this->face().i;
-        const unsigned j = this->face().j;
-        Scalar lambdaI, lambdaJ;
-
-        if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
-        {
-            lambdaI =
-                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].saturation(wPhaseIdx),
-                                                                   elemVolVars[i].thermalConductivity(wPhaseIdx),
-                                                                   elemVolVars[i].thermalConductivity(nPhaseIdx),
-                                                                   problem.spatialParams().thermalConductivitySolid(element, this->fvGeometry_, i),
-                                                                   problem.spatialParams().porosity(element, this->fvGeometry_, i));
-            lambdaJ =
-                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].saturation(wPhaseIdx),
-                                                                   elemVolVars[j].thermalConductivity(wPhaseIdx),
-                                                                   elemVolVars[j].thermalConductivity(nPhaseIdx),
-                                                                   problem.spatialParams().thermalConductivitySolid(element, this->fvGeometry_, j),
-                                                                   problem.spatialParams().porosity(element, this->fvGeometry_, j));
-        }
-        else
-        {
-            const Element& elementI = *this->fvGeometry_.neighbors[i];
-            FVElementGeometry fvGeometryI;
-            fvGeometryI.subContVol[0].global = elementI.geometry().center();
-
-            lambdaI =
-                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].saturation(wPhaseIdx),
-                                                                   elemVolVars[i].thermalConductivity(wPhaseIdx),
-                                                                   elemVolVars[i].thermalConductivity(nPhaseIdx),
-                                                                   problem.spatialParams().thermalConductivitySolid(elementI, fvGeometryI, 0),
-                                                                   problem.spatialParams().porosity(elementI, fvGeometryI, 0));
-
-            const Element& elementJ = *this->fvGeometry_.neighbors[j];
-            FVElementGeometry fvGeometryJ;
-            fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
-
-            lambdaJ =
-                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].saturation(wPhaseIdx),
-                                                                   elemVolVars[j].thermalConductivity(wPhaseIdx),
-                                                                   elemVolVars[j].thermalConductivity(nPhaseIdx),
-                                                                   problem.spatialParams().thermalConductivitySolid(elementJ, fvGeometryJ, 0),
-                                                                   problem.spatialParams().porosity(elementJ, fvGeometryJ, 0));
-        }
-
-        // -> harmonic mean
-        lambdaEff_ = harmonicMean(lambdaI, lambdaJ);
-    }
-
-private:
-    Scalar lambdaEff_;
-    Scalar normalMatrixHeatFlux_;
-    GlobalPosition temperatureGrad_;
-    int faceIdx_;
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/implicit/2p2cni/2p2cniindices.hh b/dumux/implicit/2p2cni/2p2cniindices.hh
deleted file mode 100644
index bfa42270acb1e6dfed37badbb5bc8890edd895c6..0000000000000000000000000000000000000000
--- a/dumux/implicit/2p2cni/2p2cniindices.hh
+++ /dev/null
@@ -1,49 +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 Defines the indices used by the non-isothermal two-phase two-component model
- */
-#ifndef DUMUX_2P2CNI_INDICES_HH
-#define DUMUX_2P2CNI_INDICES_HH
-
-#include "2p2cniproperties.hh"
-#include <dumux/implicit/2p2c/2p2cindices.hh>
-
-namespace Dumux
-{
-/*!
- * \ingroup TwoPTwoCNIModel
- * \ingroup ImplicitIndices
- * \brief Indices for the non-isothermal two-phase two-component model
- *
- * \tparam formulation The formulation, either pwsn or pnsw.
- * \tparam PVOffset The first index in a primary variable vector.
- */
-template <class TypeTag, int formulation, int PVOffset>
-class TwoPTwoCNIIndices : public TwoPTwoCIndices<TypeTag, formulation, PVOffset>
-{
-public:
-    static const int temperatureIdx = PVOffset + 2; //! The index for temperature in primary variable vectors.
-    static const int energyEqIdx = PVOffset + 2; //! The index for energy in equation vectors.
-};
-
-}
-#endif
diff --git a/dumux/implicit/2p2cni/2p2cnilocalresidual.hh b/dumux/implicit/2p2cni/2p2cnilocalresidual.hh
deleted file mode 100644
index 19e7a5b3202defafce0eaa6cf770ab2d7a9d79e5..0000000000000000000000000000000000000000
--- a/dumux/implicit/2p2cni/2p2cnilocalresidual.hh
+++ /dev/null
@@ -1,172 +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 Element-wise calculation of the Jacobian matrix for problems
- *        using the non-isothermal two-phase two-component fully implicit model.
- *
- */
-#ifndef DUMUX_NEW_2P2CNI_LOCAL_RESIDUAL_HH
-#define DUMUX_NEW_2P2CNI_LOCAL_RESIDUAL_HH
-
-#include "2p2cniproperties.hh"
-#include <dumux/implicit/2p2c/2p2clocalresidual.hh>
-
-namespace Dumux
-{
-/*!
- * \ingroup TwoPTwoCNIModel
- * \ingroup ImplicitLocalResidual
- * \brief Element-wise calculation of the Jacobian matrix for problems
- *        using the two-phase two-component fully implicit model.
- */
-template<class TypeTag>
-class TwoPTwoCNILocalResidual : public TwoPTwoCLocalResidual<TypeTag>
-{
-    typedef TwoPTwoCLocalResidual<TypeTag> ParentType;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-    enum {
-        numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
-        energyEqIdx = Indices::energyEqIdx,
-        temperatureIdx = Indices::temperatureIdx,
-        wPhaseIdx = Indices::wPhaseIdx,
-        nPhaseIdx = Indices::nPhaseIdx
-    };
-
-public:
-    /*!
-     * \brief Constructor. Sets the upwind weight.
-     */
-    TwoPTwoCNILocalResidual()
-    {
-        // retrieve the upwind weight for the mass conservation equations. Use the value
-        // specified via the property system as default, and overwrite
-        // it by the run-time parameter from the Dune::ParameterTree
-        massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
-    };
-
-    /*!
-     * \brief Evaluate the amount of all conservation quantities
-     *        (e.g. phase mass) within a sub-control volume.
-     *
-     * The result should be averaged over the volume (e.g. phase mass
-     * inside a sub control volume divided by the volume)
-     *
-     *  \param storage The storage of the conservation quantity (mass or energy) within the sub-control volume
-     *  \param scvIdx The sub-control-volume index
-     *  \param usePrevSol Evaluate function with solution of current or previous time step
-     */
-    void computeStorage(PrimaryVariables &storage, const int scvIdx, bool usePrevSol) const
-    {
-        // compute the storage term for phase mass
-        ParentType::computeStorage(storage, scvIdx, usePrevSol);
-
-        // if flag usePrevSol is set, the solution from the previous
-        // time step is used, otherwise the current solution is
-        // used. The secondary variables are used accordingly.  This
-        // is required to compute the derivative of the storage term
-        // using the implicit Euler method.
-        const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
-        const VolumeVariables &volVars = elemVolVars[scvIdx];
-
-        // compute the energy storage
-        storage[energyEqIdx] =
-            volVars.porosity()*(volVars.density(wPhaseIdx) *
-                                volVars.internalEnergy(wPhaseIdx) *
-                                volVars.saturation(wPhaseIdx)
-                                +
-                                volVars.density(nPhaseIdx) *
-                                volVars.internalEnergy(nPhaseIdx) *
-                                volVars.saturation(nPhaseIdx))
-            +
-            // heat capacity is already multiplied by the density
-            // of the porous material and the porosity in the problem file
-            volVars.temperature()*volVars.heatCapacity();
-    }
-
-    /*!
-     * \brief Evaluates the advective mass flux and the heat flux
-     *        over a face of a sub-control volume and writes the result into
-     *        the flux vector.
-     *
-     * \param flux The advective flux over the sub-control-volume face for each component
-     * \param fluxVars The flux variables at the current sub-control-volume face
-     *
-     * This method is called by compute flux (base class).
-     */
-    void computeAdvectiveFlux(PrimaryVariables &flux,
-                              const FluxVariables &fluxVars) const
-    {
-        // advective mass flux
-        ParentType::computeAdvectiveFlux(flux, fluxVars);
-
-        // advective heat flux in all phases
-        flux[energyEqIdx] = 0;
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-            // vertex data of the upstream and the downstream vertices
-            const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(phaseIdx));
-            const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(phaseIdx));
-
-            flux[energyEqIdx] +=
-                fluxVars.volumeFlux(phaseIdx) * (massUpwindWeight_ * // upstream vertex
-                                                 (  up.density(phaseIdx) *
-                                                    up.enthalpy(phaseIdx))
-                                                 +
-                                                 (1-massUpwindWeight_) * // downstream vertex
-                                                 (  dn.density(phaseIdx) *
-                                                    dn.enthalpy(phaseIdx)) );
-        }
-    }
-
-    /*!
-     * \brief Adds the diffusive heat flux to the flux vector over
-     *        the face of a sub-control volume.
-     *
-     * \param flux The diffusive flux over the sub-control-volume face for each conservation quantity (mass, energy)
-     * \param fluxVars The flux variables at the current sub-control-volume face
-     *
-     * This method is called by compute flux (base class).
-     */
-    void computeDiffusiveFlux(PrimaryVariables &flux,
-                              const FluxVariables &fluxVars) const
-    {
-        // diffusive mass flux
-        ParentType::computeDiffusiveFlux(flux, fluxVars);
-
-        // diffusive heat flux
-        flux[temperatureIdx] +=
-            fluxVars.normalMatrixHeatFlux();
-    }
-
-private:
-    Scalar massUpwindWeight_;
-
-};
-
-}
-
-#endif
diff --git a/dumux/implicit/2p2cni/2p2cnimodel.hh b/dumux/implicit/2p2cni/2p2cnimodel.hh
deleted file mode 100644
index 4ecaf664e73bfcf4862adc3ac72273318d009892..0000000000000000000000000000000000000000
--- a/dumux/implicit/2p2cni/2p2cnimodel.hh
+++ /dev/null
@@ -1,103 +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 Adaption of the fully implicit scheme to the non-isothermal
- *        two-phase two-component flow model.
- */
-#ifndef DUMUX_NEW_2P2CNI_MODEL_HH
-#define DUMUX_NEW_2P2CNI_MODEL_HH
-
-#include "2p2cniproperties.hh"
-#include <dumux/implicit/2p2c/2p2cmodel.hh>
-
-namespace Dumux {
-/*!
- * \ingroup TwoPTwoCNIModel
- * \brief Adaption of the fully implicit scheme to the non-isothermal
- *        two-phase two-component flow model.
- *
- * This model implements a non-isothermal two-phase flow of two compressible and partly miscible fluids
- * \f$\alpha \in \{ w, n \}\f$. Thus each component \f$\kappa \in \{ w, a \}\f$ can be present in
- * each phase.
- * Using the standard multiphase Darcy approach a mass balance equation is
- * solved:
- * \f{eqnarray*}
- && \phi \frac{\partial (\sum_\alpha \varrho_\alpha X_\alpha^\kappa S_\alpha )}{\partial t}
- - \sum_\alpha \text{div} \left\{ \varrho_\alpha X_\alpha^\kappa
- \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K}
- (\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g}) \right\}\\
- &-& \sum_\alpha \text{div} \left\{ D_{\alpha,\text{pm}}^\kappa \varrho_{\alpha} \frac{M^\kappa}{M_\alpha}
- \textbf{grad} x^\kappa_{\alpha} \right\}
- - \sum_\alpha q_\alpha^\kappa = 0 \qquad \kappa \in \{w, a\} \, ,
- \alpha \in \{w, n\}
- *     \f}
- * For the energy balance, local thermal equilibrium is assumed which results in one
- * energy conservation equation for the porous solid matrix and the fluids:
- * \f{eqnarray*}
- && \phi \frac{\partial \left( \sum_\alpha \varrho_\alpha u_\alpha S_\alpha \right)}{\partial t}
- + \left( 1 - \phi \right) \frac{\partial (\varrho_s c_s T)}{\partial t}
- - \sum_\alpha \text{div} \left\{ \varrho_\alpha h_\alpha
- \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K} \left( \textbf{grad}\,
- p_\alpha
- - \varrho_\alpha \mathbf{g} \right) \right\} \\
- &-& \text{div} \left( \lambda_\text{pm} \textbf{grad} \, T \right)
- - q^h = 0 \qquad \alpha \in \{w, n\}
- \f}
- *
- * All equations are discretized using a vertex-centered finite volume (box)
- * or cell-centered finite volume scheme as spatial
- * and the implicit Euler method as time discretization.
- *
- * By using constitutive relations for the capillary pressure \f$p_c =
- * p_n - p_w\f$ and relative permeability \f$k_{r\alpha}\f$ and taking
- * advantage of the fact that \f$S_w + S_n = 1\f$ and \f$X^\kappa_w + X^\kappa_n = 1\f$, the number of
- * unknowns can be reduced to two.
- * If both phases are present the primary variables are, like in the nonisothermal two-phase model, either \f$p_w\f$, \f$S_n\f$ and
- * temperature or \f$p_n\f$, \f$S_w\f$ and temperature. The formulation which ought to be used can be
- * specified by setting the <tt>Formulation</tt> property to either
- * <tt>TwoPTwoCIndices::pWsN</tt> or <tt>TwoPTwoCIndices::pNsW</tt>. By
- * default, the model uses \f$p_w\f$ and \f$S_n\f$.
- * In case that only one phase (nonwetting or wetting phase) is present the second primary
- * variable represents a mass fraction. The correct assignment of the second
- * primary variable is performed by a phase state dependent primary variable switch.
- * The phase state is stored for all nodes of the system. The following cases can be distinguished:
- * <ul>
- *  <li>
- *    Both phases are present: The saturation is used (either\f$S_n\f$ or \f$S_w\f$, dependent on the chosen formulation).
- *  </li>
- *  <li>
- *    Only wetting phase is present: The mass fraction of air in the wetting phase \f$X^a_w\f$ is used.
- *  </li>
- *  <li>
- *    Only non-wetting phase is present: The mass fraction of water in the non-wetting phase, \f$X^w_n\f$, is used.
- *  </li>
- * </ul>
- */
-template<class TypeTag>
-class TwoPTwoCNIModel : public TwoPTwoCModel<TypeTag>
-{
-};
-
-}
-
-#include "2p2cnipropertydefaults.hh"
-
-#endif
diff --git a/dumux/implicit/2p2cni/2p2cniproperties.hh b/dumux/implicit/2p2cni/2p2cniproperties.hh
deleted file mode 100644
index fe422ce1f2f152b53e09632c464639ee7654133f..0000000000000000000000000000000000000000
--- a/dumux/implicit/2p2cni/2p2cniproperties.hh
+++ /dev/null
@@ -1,57 +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/>.   *
- *****************************************************************************/
-/*!
- * \ingroup Properties
- * \ingroup ImplicitProperties
- * \ingroup TwoPTwoCNIModel
- * \file
- *
- * \brief Defines the properties required for the non-isothermal two-phase
- *        two-component fully implicit model.
- */
-#ifndef DUMUX_2P2CNI_PROPERTIES_HH
-#define DUMUX_2P2CNI_PROPERTIES_HH
-
-#warning You are including the old 2p2cni model which will be removed after 2.6. See CHANGELOG for details.
-
-#include <dumux/implicit/2p2c/2p2cproperties.hh>
-
-namespace Dumux
-{
-
-namespace Properties
-{
-//////////////////////////////////////////////////////////////////
-// Type tags
-//////////////////////////////////////////////////////////////////
-
-//! The type tags for the implicit non-isothermal two-phase two-component problems
-NEW_TYPE_TAG(TwoPTwoCNI, INHERITS_FROM(TwoPTwoC));
-NEW_TYPE_TAG(BoxTwoPTwoCNI, INHERITS_FROM(BoxModel, TwoPTwoCNI));
-NEW_TYPE_TAG(CCTwoPTwoCNI, INHERITS_FROM(CCModel, TwoPTwoCNI));
-
-//////////////////////////////////////////////////////////////////
-// Property tags
-//////////////////////////////////////////////////////////////////
-
-NEW_PROP_TAG(ThermalConductivityModel);   //!< The model for the effective thermal conductivity
-}
-}
-
-#endif
diff --git a/dumux/implicit/2p2cni/2p2cnipropertydefaults.hh b/dumux/implicit/2p2cni/2p2cnipropertydefaults.hh
deleted file mode 100644
index 98fedf0d485ab5eb38311298b715782d32b5ec12..0000000000000000000000000000000000000000
--- a/dumux/implicit/2p2cni/2p2cnipropertydefaults.hh
+++ /dev/null
@@ -1,87 +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/>.   *
- *****************************************************************************/
-/*!
- * \ingroup Properties
- * \ingroup ImplicitProperties
- * \ingroup TwoPTwoCNIModel
- * \file
- *
- * \brief Defines default values for most properties required by the
- *        non-isothermal two-phase two-component fully implicit model.
- */
-#ifndef DUMUX_2P2CNI_PROPERTY_DEFAULTS_HH
-#define DUMUX_2P2CNI_PROPERTY_DEFAULTS_HH
-
-#include "2p2cniproperties.hh"
-#include <dumux/implicit/2p2c/2p2cpropertydefaults.hh>
-
-#include "2p2cnimodel.hh"
-#include "2p2cniindices.hh"
-#include "2p2cnilocalresidual.hh"
-#include "2p2cnivolumevariables.hh"
-#include "2p2cnifluxvariables.hh"
-
-#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
-
-namespace Dumux
-{
-
-namespace Properties
-{
-//////////////////////////////////////////////////////////////////
-// Property values
-//////////////////////////////////////////////////////////////////
-//! Set the number of equations to 3
-SET_INT_PROP(TwoPTwoCNI, NumEq, 3); //!< set the number of equations to 3
-
-//! Use the 2p2cni local jacobian operator for the 2p2cni model
-SET_TYPE_PROP(TwoPTwoCNI,
-              LocalResidual,
-              TwoPTwoCNILocalResidual<TypeTag>);
-
-//! the Model property
-SET_TYPE_PROP(TwoPTwoCNI, Model, TwoPTwoCNIModel<TypeTag>);
-
-//! the VolumeVariables property
-SET_TYPE_PROP(TwoPTwoCNI, VolumeVariables, TwoPTwoCNIVolumeVariables<TypeTag>);
-
-
-//! the FluxVariables property
-SET_TYPE_PROP(TwoPTwoCNI, FluxVariables, TwoPTwoCNIFluxVariables<TypeTag>);
-
-//! The indices required by the non-isothermal 2p2c model
-SET_PROP(TwoPTwoCNI, Indices) 
-{ private:
-    enum { formulation = GET_PROP_VALUE(TypeTag, Formulation) };
- public:
-    typedef TwoPTwoCNIIndices<TypeTag, formulation, 0> type;
-};
-
-//! Somerton is used as default model to compute the effective thermal heat conductivity
-SET_PROP(TwoPTwoCNI, ThermalConductivityModel)
-{ private :
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-  public:
-    typedef ThermalConductivitySomerton<Scalar> type;
-};
-
-}
-
-}
-#endif
diff --git a/dumux/implicit/2p2cni/2p2cnivolumevariables.hh b/dumux/implicit/2p2cni/2p2cnivolumevariables.hh
deleted file mode 100644
index 72740d689c1eb8bfd5deb7b37e984912ab321e8f..0000000000000000000000000000000000000000
--- a/dumux/implicit/2p2cni/2p2cnivolumevariables.hh
+++ /dev/null
@@ -1,138 +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 Contains the quantities which are constant within a
- *        finite volume in the non-isothermal two-phase two-component
- *        model.
- */
-#ifndef DUMUX_2P2CNI_VOLUME_VARIABLES_HH
-#define DUMUX_2P2CNI_VOLUME_VARIABLES_HH
-
-#include "2p2cniproperties.hh"
-#include <dumux/implicit/2p2c/2p2cvolumevariables.hh>
-
-namespace Dumux
-{
-
-/*!
- * \ingroup TwoPTwoCNIModel
- * \ingroup ImplicitVolumeVariables
- * \brief Contains the quantities which are are constant within a
- *        finite volume in the non-isothermal two-phase two-component
- *        model.
- */
-template <class TypeTag>
-class TwoPTwoCNIVolumeVariables : public TwoPTwoCVolumeVariables<TypeTag>
-{
-    //! \cond false
-    typedef TwoPTwoCVolumeVariables<TypeTag> ParentType;
-    typedef typename ParentType::FluidState FluidState;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-    enum { temperatureIdx = Indices::temperatureIdx };
-
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-    typedef typename GridView::template Codim<0>::Entity Element;
-    //! \endcond
-
-public:
-    /*!
-     * \brief Returns the total internal energy of a phase in the
-     *        sub-control volume in \f$\mathrm{[J/kg]}\f$.
-     *
-     * \param phaseIdx The phase index
-     */
-    Scalar internalEnergy(const int phaseIdx) const
-    { return this->fluidState_.internalEnergy(phaseIdx); };
-
-    /*!
-     * \brief Returns the total enthalpy of a phase in the sub-control
-     *        volume in \f$\mathrm{[J/kg]}\f$.
-     *
-     * \param phaseIdx The phase index
-     */
-    Scalar enthalpy(const int phaseIdx) const
-    { return this->fluidState_.enthalpy(phaseIdx); };
-
-    /*!
-     * \brief Returns the total heat capacity \f$\mathrm{[J/(K*m^3]}\f$ of the rock matrix in
-     *        the sub-control volume.
-     */
-    Scalar heatCapacity() const
-    { return heatCapacity_; };
-
-    /*!
-     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ of the fluid phase in
-     *        the sub-control volume.
-     */
-    Scalar thermalConductivity(const int phaseIdx) const
-    { return FluidSystem::thermalConductivity(this->fluidState_, phaseIdx); };
-
-
-protected:
-    // this method gets called by the parent class. since this method
-    // is protected, we are friends with our parent..
-    friend class TwoPTwoCVolumeVariables<TypeTag>;
-
-    static Scalar temperature_(const PrimaryVariables &priVars,
-                               const Problem& problem,
-                               const Element &element,
-                               const FVElementGeometry &fvGeometry,
-                               const int scvIdx)
-    {
-        return priVars[temperatureIdx];
-    }
-
-    template<class ParameterCache>
-    static Scalar enthalpy_(const FluidState& fluidState,
-                            const ParameterCache& paramCache,
-                            const int phaseIdx)
-    {
-        return FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
-    }
-
-    /*!
-     * \copydoc ImplicitVolumeVariables::update
-     */
-    void updateEnergy_(const PrimaryVariables &priVars,
-                       const Problem &problem,
-                       const Element &element,
-                       const FVElementGeometry &fvGeometry,
-                       const int scvIdx,
-                       bool isOldSol)
-    {
-        // compute and set the heat capacity of the solid phase
-        heatCapacity_ = problem.spatialParams().heatCapacity(element, fvGeometry, scvIdx);
-        Valgrind::CheckDefined(heatCapacity_);
-    };
-
-    Scalar heatCapacity_;
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/implicit/2p2cni/Makefile.am b/dumux/implicit/2p2cni/Makefile.am
deleted file mode 100644
index 8fa2ae7c42315b387457e4059305b851f949b6c0..0000000000000000000000000000000000000000
--- a/dumux/implicit/2p2cni/Makefile.am
+++ /dev/null
@@ -1,4 +0,0 @@
-2p2cnidir = $(includedir)/dumux/implicit/2p2cni
-2p2cni_HEADERS := $(wildcard *.hh)
-
-include $(top_srcdir)/am/global-rules
diff --git a/dumux/implicit/2pni/2pnifluxvariables.hh b/dumux/implicit/2pni/2pnifluxvariables.hh
deleted file mode 100644
index e0d91745f3fbcf05159ef729a0f8719ce50f2fac..0000000000000000000000000000000000000000
--- a/dumux/implicit/2pni/2pnifluxvariables.hh
+++ /dev/null
@@ -1,208 +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 data which is required to calculate
- *        the heat fluxes over a face of a finite volume for the non-isothermal
- *        two-phase model.
- *
- * This means pressure and temperature gradients, phase densities at
- * the integration point, etc.
- */
-#ifndef DUMUX_2PNI_FLUX_VARIABLES_HH
-#define DUMUX_2PNI_FLUX_VARIABLES_HH
-
-#include "2pniproperties.hh"
-#include <dumux/common/math.hh>
-#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
-
-namespace Dumux
-{
-
-/*!
- * \ingroup TwoPNIModel
- * \ingroup ImplicitFluxVariables
- * \brief This file contains data which is required to calculate
- *        the heat fluxes over a face of a finite volume for the non-isothermal
- *        two-phase model.
- */
-template <class TypeTag>
-class TwoPNIFluxVariables : public ImplicitDarcyFluxVariables<TypeTag>
-{
-    typedef ImplicitDarcyFluxVariables<TypeTag> ParentType;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-    typedef typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel) ThermalConductivityModel;
-
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-    typedef typename GridView::template Codim<0>::Entity Element;
-    enum { dim = GridView::dimension };
-    enum { dimWorld = GridView::dimensionworld };
-
-    typedef typename GridView::ctype CoordScalar;
-    typedef Dune::FieldVector<CoordScalar, dim> DimVector;
-    typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-    enum {
-        wPhaseIdx = Indices::wPhaseIdx,
-        nPhaseIdx = Indices::nPhaseIdx
-    };
-
-public:
-
-    /*!
-     * \brief The constructor
-     *
-     * \param problem The problem
-     * \param element The finite element
-     * \param fvGeometry The finite-volume geometry in the fully implicit scheme
-     * \param faceIdx The local index of the SCV (sub-control-volume) face
-     * \param elemVolVars The volume variables of the current element
-     * \param onBoundary Distinguishes if we are on a sub-control-volume face or on a boundary face
-     */
-
-    TwoPNIFluxVariables(const Problem &problem,
-                        const Element &element,
-                        const FVElementGeometry &fvGeometry,
-                        int faceIdx,
-                        const ElementVolumeVariables &elemVolVars,
-                        const bool onBoundary = false)
-        : ParentType(problem, element, fvGeometry, faceIdx, elemVolVars, onBoundary)
-    {
-        faceIdx_ = faceIdx;
-
-        calculateValues_(problem, element, elemVolVars);
-    }
-
-    /*!
-     * \brief The total heat flux \f$\mathrm{[J/s]}\f$ due to heat conduction
-     *        of the rock matrix over the sub-control volume's face in
-     *        direction of the face normal.
-     */
-    Scalar normalMatrixHeatFlux() const
-    { return normalMatrixHeatFlux_; }
-
-
-    /*!
-     * \brief The local temperature gradient at the IP of the considered scv face.
-     */
-    GlobalPosition temperatureGradient() const
-    { return temperatureGrad_; }
-
-    /*!
-     * \brief The harmonically averaged effective thermal conductivity.
-     */
-    Scalar effThermalConductivity() const
-    { return lambdaEff_; }
-
-protected:
-    void calculateValues_(const Problem &problem,
-                          const Element &element,
-                          const ElementVolumeVariables &elemVolVars)
-    {
-        // calculate temperature gradient using finite element
-        // gradients
-        temperatureGrad_ = 0;
-        GlobalPosition tmp(0.0);
-        for (unsigned int idx = 0; idx < this->face().numFap; idx++)
-        {
-            tmp = this->face().grad[idx];
-
-            // index for the element volume variables
-            int volVarsIdx = this->face().fapIndices[idx];
-
-            tmp *= elemVolVars[volVarsIdx].temperature();
-            temperatureGrad_ += tmp;
-        }
-
-        lambdaEff_ = 0;
-        calculateEffThermalConductivity_(problem, element, elemVolVars);
-
-        // project the heat flux vector on the face's normal vector
-        normalMatrixHeatFlux_ = temperatureGrad_ * this->face().normal;
-        normalMatrixHeatFlux_ *= -lambdaEff_;
-    }
-
-    void calculateEffThermalConductivity_(const Problem &problem,
-                                          const Element &element,
-                                          const ElementVolumeVariables &elemVolVars)
-    {
-        const unsigned i = this->face().i;
-        const unsigned j = this->face().j;
-        Scalar lambdaI, lambdaJ;
-
-        if (GET_PROP_VALUE(TypeTag, ImplicitIsBox))
-        {
-            lambdaI =
-                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].saturation(wPhaseIdx),
-                                                                   elemVolVars[i].thermalConductivity(wPhaseIdx),
-                                                                   elemVolVars[i].thermalConductivity(nPhaseIdx),
-                                                                   problem.spatialParams().thermalConductivitySolid(element, this->fvGeometry_, i),
-                                                                   problem.spatialParams().porosity(element, this->fvGeometry_, i));
-            lambdaJ =
-                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].saturation(wPhaseIdx),
-                                                                   elemVolVars[j].thermalConductivity(wPhaseIdx),
-                                                                   elemVolVars[j].thermalConductivity(nPhaseIdx),
-                                                                   problem.spatialParams().thermalConductivitySolid(element, this->fvGeometry_, j),
-                                                                   problem.spatialParams().porosity(element, this->fvGeometry_, j));
-        }
-        else
-        {
-            const Element& elementI = *this->fvGeometry_.neighbors[i];
-            FVElementGeometry fvGeometryI;
-            fvGeometryI.subContVol[0].global = elementI.geometry().center();
-
-            lambdaI =
-                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[i].saturation(wPhaseIdx),
-                                                                   elemVolVars[i].thermalConductivity(wPhaseIdx),
-                                                                   elemVolVars[i].thermalConductivity(nPhaseIdx),
-                                                                   problem.spatialParams().thermalConductivitySolid(elementI, fvGeometryI, 0),
-                                                                   problem.spatialParams().porosity(elementI, fvGeometryI, 0));
-
-            const Element& elementJ = *this->fvGeometry_.neighbors[j];
-            FVElementGeometry fvGeometryJ;
-            fvGeometryJ.subContVol[0].global = elementJ.geometry().center();
-
-            lambdaJ =
-                ThermalConductivityModel::effectiveThermalConductivity(elemVolVars[j].saturation(wPhaseIdx),
-                                                                   elemVolVars[j].thermalConductivity(wPhaseIdx),
-                                                                   elemVolVars[j].thermalConductivity(nPhaseIdx),
-                                                                   problem.spatialParams().thermalConductivitySolid(elementJ, fvGeometryJ, 0),
-                                                                   problem.spatialParams().porosity(elementJ, fvGeometryJ, 0));
-        }
-
-        // -> harmonic mean
-        lambdaEff_ = harmonicMean(lambdaI, lambdaJ);
-    }
-
-
-private:
-    Scalar lambdaEff_;
-    Scalar normalMatrixHeatFlux_;
-    GlobalPosition temperatureGrad_;
-    int faceIdx_;
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/implicit/2pni/2pniindices.hh b/dumux/implicit/2pni/2pniindices.hh
deleted file mode 100644
index c0e48141b6ff54b9d142c7445ce73e096df5e00a..0000000000000000000000000000000000000000
--- a/dumux/implicit/2pni/2pniindices.hh
+++ /dev/null
@@ -1,50 +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 Defines the indices used by the non-isothermal two-phase fully implicit model.
- */
-#ifndef DUMUX_2PNI_INDICES_HH
-#define DUMUX_2PNI_INDICES_HH
-
-#include "2pniproperties.hh"
-#include <dumux/implicit/2p/2pindices.hh>
-
-namespace Dumux
-{
-// \{
-
-/*!
- * \ingroup TwoPNIModel
- * \ingroup ImplicitIndices
- * \brief Defines the indices used by the non-isothermal two-phase fully implicit model.
- *
- * \tparam PVOffset The first index in a primary variable vector.
- */
-template <class TypeTag, int PVOffset = 0>
-class TwoPNIIndices : public TwoPIndices<TypeTag, PVOffset>
-{
-public:
-    static const int temperatureIdx = PVOffset + 2; //! The primary variable index for temperature
-    static const int energyEqIdx = PVOffset + 2; //! The equation index of the energy equation
-};
-// \}
-}
-
-#endif
diff --git a/dumux/implicit/2pni/2pnilocalresidual.hh b/dumux/implicit/2pni/2pnilocalresidual.hh
deleted file mode 100644
index 705fbd1caa2b2371260cc65384c68eb9785e03e8..0000000000000000000000000000000000000000
--- a/dumux/implicit/2pni/2pnilocalresidual.hh
+++ /dev/null
@@ -1,175 +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 Element-wise calculation of the Jacobian matrix for problems
- *        using the non-isothermal two-phase fully implicit model.
- *
- */
-#ifndef DUMUX_NEW_2PNI_LOCAL_RESIDUAL_HH
-#define DUMUX_NEW_2PNI_LOCAL_RESIDUAL_HH
-
-#include "2pniproperties.hh"
-#include <dumux/implicit/2p/2plocalresidual.hh>
-
-namespace Dumux
-{
-/*!
- * \ingroup TwoPNIModel
- * \ingroup ImplicitLocalResidual
- * \brief Element-wise calculation of the Jacobian matrix for problems
- *        using the non-isothermal two-phase fully implicit model.
- */
-
-template<class TypeTag>
-class TwoPNILocalResidual : public TwoPLocalResidual<TypeTag>
-{
-    typedef TwoPLocalResidual<TypeTag> ParentType;
-
-    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-    enum {
-        numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
-        temperatureIdx = Indices::temperatureIdx,
-        energyEqIdx = Indices::energyEqIdx,
-        wPhaseIdx = Indices::wPhaseIdx,
-        nPhaseIdx = Indices::nPhaseIdx
-    };
-
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-
-public:
-    /*!
-     * \brief Constructor. Sets the upwind weight.
-     */
-    TwoPNILocalResidual()
-    {
-        // retrieve the upwind weight for the mass conservation equations. Use the value
-        // specified via the property system as default, and overwrite
-        // it by the run-time parameter from the Dune::ParameterTree
-        massUpwindWeight_ = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
-    };
-
-    /*!
-     * \brief Evaluate the amount all conservation quantites
-     *        (e.g. phase mass and energy storage) within a sub-control volume.
-     *
-     * The result should be averaged over the volume (e.g. phase mass
-     * inside a sub control volume divided by the volume)
-     *
-     *  \param storage The phase mass within the sub-control volume
-     *  \param scvIdx The sub-control volume index
-     *  \param usePrevSol Evaluate function with solution of current or previous time step
-     */
-    void computeStorage(PrimaryVariables &storage, int scvIdx, bool usePrevSol) const
-    {
-        // compute the storage term for phase mass
-        ParentType::computeStorage(storage, scvIdx, usePrevSol);
-
-        // if flag usePrevSol is set, the solution from the previous
-        // time step is used, otherwise the current solution is
-        // used. The secondary variables are used accordingly.  This
-        // is required to compute the derivative of the storage term
-        // using the implicit euler method.
-        const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
-        const VolumeVariables &volVars = elemVolVars[scvIdx];
-
-        // compute the energy storage
-        storage[temperatureIdx] =
-            volVars.porosity()*(volVars.density(wPhaseIdx) *
-                                volVars.internalEnergy(wPhaseIdx) *
-                                volVars.saturation(wPhaseIdx)
-                                +
-                                volVars.density(nPhaseIdx) *
-                                volVars.internalEnergy(nPhaseIdx) *
-                                volVars.saturation(nPhaseIdx))
-          + volVars.temperature()*volVars.heatCapacity();
-    }
-
-    /*!
-     * \brief Evaluates the advective mass flux and the heat flux
-     * over a face of a subcontrol volume and writes the result in
-     * the flux vector.
-     *
-     * \param flux The advective flux over the sub-control-volume face for each phase
-     * \param fluxVars The flux variables at the current SCV
-     *
-     *
-     * This method is called by compute flux (base class)
-     */
-    void computeAdvectiveFlux(PrimaryVariables &flux,
-                              const FluxVariables &fluxVars) const
-    {
-        // advective mass flux
-        ParentType::computeAdvectiveFlux(flux, fluxVars);
-
-        // advective heat flux in all phases
-        flux[energyEqIdx] = 0;
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-            // data attached to upstream and the downstream vertices
-            // of the current phase
-            const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(phaseIdx));
-            const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(phaseIdx));
-
-            // add advective energy flux in current phase
-            flux[energyEqIdx] +=
-                fluxVars.volumeFlux(phaseIdx)
-                *
-                ((    massUpwindWeight_)*
-                 up.density(phaseIdx)*
-                 up.enthalpy(phaseIdx)
-                 +
-                 (1 - massUpwindWeight_)*
-                 dn.density(phaseIdx)*
-                 dn.enthalpy(phaseIdx));
-        }
-    }
-
-    /*!
-     * \brief Adds the diffusive heat flux to the flux vector over
-     *        the face of a sub-control volume.
-     *
-     * \param flux The diffusive flux over the sub-control-volume face for each phase
-     * \param fluxVars The flux variables at the current SCV
-     *
-     * This method is called by compute flux (base class).
-     */
-    void computeDiffusiveFlux(PrimaryVariables &flux,
-                              const FluxVariables &fluxVars) const
-    {
-        // diffusive mass flux
-        ParentType::computeDiffusiveFlux(flux, fluxVars);
-
-        // diffusive heat flux
-        flux[energyEqIdx] += fluxVars.normalMatrixHeatFlux();
-    }
-
-private:
-    Scalar massUpwindWeight_;
-
-};
-
-}
-
-#endif
diff --git a/dumux/implicit/2pni/2pnimodel.hh b/dumux/implicit/2pni/2pnimodel.hh
deleted file mode 100644
index 8be48fa50debed274091500e684f6927df603341..0000000000000000000000000000000000000000
--- a/dumux/implicit/2pni/2pnimodel.hh
+++ /dev/null
@@ -1,94 +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 Adaption of the fully implicit scheme to the non-isothermal two-phase flow model.
- */
-#ifndef DUMUX_2PNI_MODEL_HH
-#define DUMUX_2PNI_MODEL_HH
-
-#include "2pniproperties.hh"
-#include <dumux/implicit/2p/2pmodel.hh>
-
-namespace Dumux {
-
-/*!
- * \ingroup TwoPNIModel
- * \brief Adaption of the fully implicit scheme to the non-isothermal two-phase flow model.
- *
- * This model implements a non-isothermal two-phase flow for two
- * immiscible fluids \f$\alpha \in \{ w, n \}\f$. Using the standard
- * multiphase Darcy approach, the mass conservation equations for both
- * phases can be described as follows:
- * \f[
- \phi \frac{\partial \varrho_\alpha S_\alpha}{\partial t}
- - 
- \text{div} 
- \left\{ 
- \varrho_\alpha \frac{k_{r\alpha}}{\mu_\alpha} \textbf{K}
- \left( \textbf{grad}\, p_\alpha - \varrho_{\alpha} \mathbf{g} \right)
- \right\}
- -
- q_\alpha = 0 \qquad \alpha \in \{w, n\}
- \f]
- *
- * For the energy balance, local thermal equilibrium is assumed. This
- * results in one energy conservation equation for the porous solid
- * matrix and the fluids: 
- 
- \f{align*}{
- \phi \frac{\partial \sum_\alpha \varrho_\alpha u_\alpha S_\alpha}{\partial t}
- & + 
- \left( 1 - \phi \right) \frac{\partial (\varrho_s c_s T)}{\partial t}
- - 
- \sum_\alpha \text{div} 
- \left\{
- \varrho_\alpha h_\alpha
- \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K} 
- \left( \textbf{grad}\,p_\alpha - \varrho_\alpha \mbox{\bf g} \right)
- \right\} \\
-    & - \text{div} \left(\lambda_{pm} \textbf{grad} \, T \right)
-    - q^h = 0, \qquad \alpha \in \{w, n\} \;,
- \f}
- * where \f$h_\alpha\f$ is the specific enthalpy of a fluid phase
- * \f$\alpha\f$ and \f$u_\alpha = h_\alpha -
- * p_\alpha/\varrho_\alpha\f$ is the specific internal energy of the
- * phase.
- *
- * All equations are discretized using a vertex-centered finite volume (box)
- * or cell-centered finite volume scheme as spatial
- * and the implicit Euler method as time discretization.
- *
- * Currently the model supports choosing either \f$p_w\f$, \f$S_n\f$
- * and \f$T\f$ or \f$p_n\f$, \f$S_w\f$ and \f$T\f$ as primary
- * variables. The formulation which ought to be used can be specified
- * by setting the <tt>Formulation</tt> property to either
- * <tt>TwoPNIIndices::pWsN</tt> or <tt>TwoPIndices::pNsW</tt>. By
- * default, the model uses \f$p_w\f$, \f$S_n\f$ and \f$T\f$.
- */
-template<class TypeTag>
-class TwoPNIModel: public TwoPModel<TypeTag>
-{};
-
-}
-
-#include "2pnipropertydefaults.hh"
-
-#endif
diff --git a/dumux/implicit/2pni/2pniproperties.hh b/dumux/implicit/2pni/2pniproperties.hh
deleted file mode 100644
index cc570df0af32f629e835f9e4f11688e683913e99..0000000000000000000000000000000000000000
--- a/dumux/implicit/2pni/2pniproperties.hh
+++ /dev/null
@@ -1,57 +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/>.   *
- *****************************************************************************/
-/*!
- * \ingroup Properties
- * \ingroup ImplicitProperties
- *  \ingroup TwoPNIModel
- * \file
- *
- * \brief Defines the properties required for the non-isotherm two-phase fully implicit model.
- */
-#ifndef DUMUX_2PNI_PROPERTIES_HH
-#define DUMUX_2PNI_PROPERTIES_HH
-
-#warning You are including the old 2pni model which will be removed after 2.6. See CHANGELOG for details.
-
-#include <dumux/implicit/2p/2pproperties.hh>
-
-namespace Dumux
-{
-
-namespace Properties
-{
-//////////////////////////////////////////////////////////////////
-// Type tags
-//////////////////////////////////////////////////////////////////
-
-//! The type tags for the implicit non-isothermal two-phase problems
-NEW_TYPE_TAG(TwoPNI, INHERITS_FROM(TwoP));
-NEW_TYPE_TAG(BoxTwoPNI, INHERITS_FROM(BoxModel, TwoPNI));
-NEW_TYPE_TAG(CCTwoPNI, INHERITS_FROM(CCModel, TwoPNI));
-
-//////////////////////////////////////////////////////////////////
-// Property tags
-//////////////////////////////////////////////////////////////////
-
-NEW_PROP_TAG(ThermalConductivityModel);   //!< The model for the effective thermal conductivity
-
-}
-}
-
-#endif
diff --git a/dumux/implicit/2pni/2pnipropertydefaults.hh b/dumux/implicit/2pni/2pnipropertydefaults.hh
deleted file mode 100644
index ddcdac848ccf8f5833de9aeb98311a63943530fb..0000000000000000000000000000000000000000
--- a/dumux/implicit/2pni/2pnipropertydefaults.hh
+++ /dev/null
@@ -1,82 +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/>.   *
- *****************************************************************************/
-/*!
- * \ingroup Properties
- * \ingroup ImplicitProperties
- * \ingroup TwoPNIModel
- * \file
- *
- * \brief Defines the default values for most of the properties
- *        required by the non-isothermal two-phase fully implicit model.
- */
-
-#ifndef DUMUX_2PNI_PROPERTY_DEFAULTS_HH
-#define DUMUX_2PNI_PROPERTY_DEFAULTS_HH
-
-#include "2pniproperties.hh"
-#include "2pnimodel.hh"
-#include "2pnilocalresidual.hh"
-#include "2pnivolumevariables.hh"
-#include "2pnifluxvariables.hh"
-#include "2pniindices.hh"
-
-#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
-
-namespace Dumux
-{
-
-
-namespace Properties
-{
-//////////////////////////////////////////////////////////////////
-// Property values
-//////////////////////////////////////////////////////////////////
-//! Set the number of equations to 3
-SET_INT_PROP(TwoPNI, NumEq, 3);
-
-//! Use the 2pni local jacobian operator for the 2pni model
-SET_TYPE_PROP(TwoPNI,
-              LocalResidual,
-              TwoPNILocalResidual<TypeTag>);
-
-//! the Model property
-SET_TYPE_PROP(TwoPNI, Model, TwoPNIModel<TypeTag>);
-
-//! the VolumeVariables property
-SET_TYPE_PROP(TwoPNI, VolumeVariables, TwoPNIVolumeVariables<TypeTag>);
-
-//! the FluxVariables property
-SET_TYPE_PROP(TwoPNI, FluxVariables, TwoPNIFluxVariables<TypeTag>);
-
-//! The indices required by the non-isothermal two-phase model
-SET_TYPE_PROP(TwoPNI, Indices, TwoPNIIndices<TypeTag, 0>);
-
-//! Somerton is used as default model to compute the effective thermal heat conductivity
-SET_PROP(TwoPNI, ThermalConductivityModel)
-{ private :
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-  public:
-    typedef ThermalConductivitySomerton<Scalar> type;
-};
-
-}
-
-}
-
-#endif
diff --git a/dumux/implicit/2pni/2pnivolumevariables.hh b/dumux/implicit/2pni/2pnivolumevariables.hh
deleted file mode 100644
index 9c3fab9dc7791d1c91223d5c4af605872d511598..0000000000000000000000000000000000000000
--- a/dumux/implicit/2pni/2pnivolumevariables.hh
+++ /dev/null
@@ -1,133 +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 Contains the quantities which are are constant within a
- *        finite volume in the non-isothermal two-phase model.
- */
-#ifndef DUMUX_2PNI_VOLUME_VARIABLES_HH
-#define DUMUX_2PNI_VOLUME_VARIABLES_HH
-
-#include "2pniproperties.hh"
-#include <dumux/implicit/2p/2pvolumevariables.hh>
-
-namespace Dumux
-{
-
-/*!
- * \ingroup TwoPNIModel
- * \ingroup ImplicitVolumeVariables
- * \brief Contains the quantities which are constant within a
- *        finite volume in the non-isothermal two-phase model.
- */
-template <class TypeTag>
-class TwoPNIVolumeVariables : public TwoPVolumeVariables<TypeTag>
-{
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
-    typedef typename GET_PROP_TYPE(TypeTag, FluidState) FluidState;
-    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-    enum { temperatureIdx = Indices::temperatureIdx };
-
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-    typedef typename GridView::template Codim<0>::Entity Element;
-
-public:
-    /*!
-     * \brief Returns the total internal energy of a phase in the
-     *        sub-control volume in \f$\mathrm{[J/kg]}\f$.
-     *
-     * \param phaseIdx The phase index
-     *
-     */
-    Scalar internalEnergy(int phaseIdx) const
-    { return this->fluidState_.internalEnergy(phaseIdx); };
-
-    /*!
-     * \brief Returns the total enthalpy of a phase in the sub-control
-     *        volume in \f$\mathrm{[J/kg]}\f$.
-     *
-     *  \param phaseIdx The phase index
-     */
-    Scalar enthalpy(int phaseIdx) const
-    { return this->fluidState_.enthalpy(phaseIdx); };
-
-    /*!
-     * \brief Returns the total heat capacity \f$\mathrm{[J/K*m^3]}\f$ of the rock matrix in
-     *        the sub-control volume.
-     */
-    Scalar heatCapacity() const
-    { return heatCapacity_; };
-
-    /*!
-     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ of the fluid phase in
-     *        the sub-control volume.
-     */
-    Scalar thermalConductivity(const int phaseIdx) const
-    { return FluidSystem::thermalConductivity(this->fluidState_, phaseIdx); };
-
-
-protected:
-    // this method gets called by the parent class. since this method
-    // is protected, we are friends with our parent..
-    friend class TwoPVolumeVariables<TypeTag>;
-
-    static Scalar temperature_(const PrimaryVariables &priVars,
-                            const Problem& problem,
-                            const Element &element,
-                            const FVElementGeometry &fvGeometry,
-                            int scvIdx)
-    {
-        return priVars[Indices::temperatureIdx];
-    }
-
-    template<class ParameterCache>
-    static Scalar enthalpy_(const FluidState& fluidState,
-                            const ParameterCache& paramCache,
-                            int phaseIdx)
-    {
-        return FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
-    }
-
-    /*!
-     * \brief Called by update() to compute the energy related quantities
-     */
-    void updateEnergy_(const PrimaryVariables &priVars,
-                       const Problem &problem,
-                       const Element &element,
-                       const FVElementGeometry &fvGeometry,
-                       int scvIdx,
-                       bool isOldSol)
-    {
-        // copmute and set the heat capacity of the solid phase
-        heatCapacity_ = problem.spatialParams().heatCapacity(element, fvGeometry, scvIdx);
-        Valgrind::CheckDefined(heatCapacity_);
-    }
-
-    Scalar heatCapacity_;
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/implicit/2pni/Makefile.am b/dumux/implicit/2pni/Makefile.am
deleted file mode 100644
index 62ed49a1405cbc952a66b654d64d78cb6ee889c5..0000000000000000000000000000000000000000
--- a/dumux/implicit/2pni/Makefile.am
+++ /dev/null
@@ -1,4 +0,0 @@
-2pnidir = $(includedir)/dumux/implicit/2pni
-2pni_HEADERS := $(wildcard *.hh)
-
-include $(top_srcdir)/am/global-rules
diff --git a/dumux/implicit/3p3c/3p3cproperties.hh b/dumux/implicit/3p3c/3p3cproperties.hh
index c2b8e0e4a48bf3a7937ce852bef60d994eae4022..81a97b4c97a333884be0bc8bbaa75f96ca88325e 100644
--- a/dumux/implicit/3p3c/3p3cproperties.hh
+++ b/dumux/implicit/3p3c/3p3cproperties.hh
@@ -49,12 +49,10 @@ NEW_TYPE_TAG(ThreePThreeC);
 NEW_TYPE_TAG(BoxThreePThreeC, INHERITS_FROM(BoxModel, ThreePThreeC));
 NEW_TYPE_TAG(CCThreePThreeC, INHERITS_FROM(CCModel, ThreePThreeC));
 
-#ifndef DUMUX_3P3CNI_PROPERTIES_HH
 //! The type tags for the corresponding non-isothermal problems
 NEW_TYPE_TAG(ThreePThreeCNI, INHERITS_FROM(ThreePThreeC, NonIsothermal));
 NEW_TYPE_TAG(BoxThreePThreeCNI, INHERITS_FROM(BoxModel, ThreePThreeCNI));
 NEW_TYPE_TAG(CCThreePThreeCNI, INHERITS_FROM(CCModel, ThreePThreeCNI));
-#endif
 
 //////////////////////////////////////////////////////////////////
 // Property tags
diff --git a/dumux/implicit/3p3cni/3p3cnifluxvariables.hh b/dumux/implicit/3p3cni/3p3cnifluxvariables.hh
deleted file mode 100644
index 541025cc60811831d7181b74995b6d2a91f26879..0000000000000000000000000000000000000000
--- a/dumux/implicit/3p3cni/3p3cnifluxvariables.hh
+++ /dev/null
@@ -1,130 +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 data which is required to calculate
- *        the heat fluxes over a face of a finite volume.
- *
- * This means temperature gradients and the normal matrix
- * heat flux.
- */
-#ifndef DUMUX_3P3CNI_FLUX_VARIABLES_HH
-#define DUMUX_3P3CNI_FLUX_VARIABLES_HH
-
-#include <dumux/common/math.hh>
-#include <dumux/implicit/3p3c/3p3cfluxvariables.hh>
-
-namespace Dumux
-{
-
-/*!
- * \ingroup ThreePThreeCNIModel
- * \ingroup ImplicitFluxVariables
- * \brief This template class contains data which is required to
- *        calculate the heat fluxes over a face of a finite
- *        volume for the non-isothermal three-phase three-component model.
- *        The mass fluxes are computed in the parent class.
- *
- * This means temperature gradients and the normal matrix
- * heat flux.
- */
-template <class TypeTag>
-class ThreePThreeCNIFluxVariables : public ThreePThreeCFluxVariables<TypeTag>
-{
-    typedef ThreePThreeCFluxVariables<TypeTag> ParentType;
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-
-    typedef typename GridView::ctype CoordScalar;
-    typedef typename GridView::template Codim<0>::Entity Element;
-    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
-
-    enum {
-        dim = GridView::dimension,
-        dimWorld = GridView::dimensionworld
-    };
-
-    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-
-    typedef Dune::FieldVector<CoordScalar, dim> DimVector;
-    typedef Dune::FieldVector<CoordScalar, dimWorld> GlobalPosition;
-
-public:
-    /*!
-     * \brief The constructor
-     *
-     * \param problem The problem
-     * \param element The finite element
-     * \param fvGeometry The finite-volume geometry in the fully implicit scheme
-     * \param faceIdx The local index of the SCV (sub-control-volume) face
-     * \param elemVolVars The volume variables of the current element
-     * \param onBoundary Distinguishes if we are on a sub-control-volume face or on a boundary face
-     */
-    ThreePThreeCNIFluxVariables(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)
-    {
-        // calculate temperature gradient using finite element
-        // gradients
-        GlobalPosition temperatureGrad(0);
-        GlobalPosition tmp(0.0);
-        for (unsigned int idx = 0; idx < this->face().numFap; idx++)
-        {
-            tmp = this->face().grad[idx];
-
-            // index for the element volume variables 
-            int volVarsIdx = this->face().fapIndices[idx];
-
-            tmp *= elemVolVars[volVarsIdx].temperature();
-            temperatureGrad += tmp;
-        }
-
-        // The spatial parameters calculates the actual heat flux vector
-        problem.spatialParams().matrixHeatFlux(tmp,
-                                                   *this,
-                                                   elemVolVars,
-                                                   temperatureGrad,
-                                                   element,
-                                                   fvGeometry,
-                                                   faceIdx);
-        // project the heat flux vector on the face's normal vector
-        normalMatrixHeatFlux_ = tmp*this->face().normal;
-    }
-
-    /*!
-     * \brief The total heat flux \f$\mathrm{[J/s]}\f$ due to heat conduction
-     *        of the rock matrix over the sub-control volume's face in
-     *        direction of the face normal.
-     */
-    Scalar normalMatrixHeatFlux() const
-    { return normalMatrixHeatFlux_; }
-
-private:
-    Scalar normalMatrixHeatFlux_;
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/implicit/3p3cni/3p3cniindices.hh b/dumux/implicit/3p3cni/3p3cniindices.hh
deleted file mode 100644
index 8dc189691ccdd4ff677007fa36c04460c1a08441..0000000000000000000000000000000000000000
--- a/dumux/implicit/3p3cni/3p3cniindices.hh
+++ /dev/null
@@ -1,47 +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 Defines the indices used by the non-isothermal three-phase three-component model
- */
-#ifndef DUMUX_3P3CNI_INDICES_HH
-#define DUMUX_3P3CNI_INDICES_HH
-
-#include <dumux/implicit/3p3c/3p3cindices.hh>
-
-namespace Dumux
-{
-/*!
- * \ingroup ThreePThreeCNIModel
- * \ingroup ImplicitIndices
- * \brief Enumerations for the non-isothermal three-phase three-component model
- *
- * \tparam PVOffset The first index in a primary variable vector.
- */
-template <class TypeTag, int PVOffset>
-class ThreePThreeCNIIndices : public ThreePThreeCIndices<TypeTag, PVOffset>
-{
-public:
-    static const int temperatureIdx = PVOffset + 3; //! The index for temperature in primary variable vectors.
-    static const int energyEqIdx = PVOffset + 3; //! The index for energy in equation vectors.
-};
-
-}
-#endif
diff --git a/dumux/implicit/3p3cni/3p3cnilocalresidual.hh b/dumux/implicit/3p3cni/3p3cnilocalresidual.hh
deleted file mode 100644
index 518304918d25b81cbad6c5b63925cc8e9a9bfbc8..0000000000000000000000000000000000000000
--- a/dumux/implicit/3p3cni/3p3cnilocalresidual.hh
+++ /dev/null
@@ -1,172 +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 Element-wise calculation of the residual for problems
- *        using the non-isothermal three-phase three-component fully implicit model.
- *
- */
-#ifndef DUMUX_NEW_3P3CNI_LOCAL_RESIDUAL_HH
-#define DUMUX_NEW_3P3CNI_LOCAL_RESIDUAL_HH
-
-#include "3p3cniproperties.hh"
-#include <dumux/implicit/3p3c/3p3clocalresidual.hh>
-
-namespace Dumux
-{
-/*!
- * \ingroup ThreePThreeCNIModel
- * \ingroup ImplicitLocalResidual
- * \brief Element-wise calculation of the residual for problems
- *        using the non-isothermal three-phase three-component fully implicit model.
- */
-template<class TypeTag>
-class ThreePThreeCNILocalResidual : public ThreePThreeCLocalResidual<TypeTag>
-{
-    typedef ThreePThreeCLocalResidual<TypeTag> ParentType;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-
-
-    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-
-    enum {
-
-        numPhases = GET_PROP_VALUE(TypeTag, NumPhases),
-
-        energyEqIdx = Indices::energyEqIdx,
-        temperatureIdx = Indices::temperatureIdx,
-
-        wPhaseIdx = Indices::wPhaseIdx,
-        nPhaseIdx = Indices::nPhaseIdx,
-        gPhaseIdx = Indices::gPhaseIdx
-    };
-
-
-    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
-    
-public:
-    /*!
-     * \brief Evaluate the amount of all conservation quantities
-     *        (e.g. phase mass) within a sub-control volume.
-     *
-     * The result should be averaged over the volume (e.g. phase mass
-     * inside a sub control volume divided by the volume)
-     *
-     *  \param storage The storage of the conservation quantitiy (mass or energy) within the sub-control volume
-     *  \param scvIdx The sub-control-volume index
-     *  \param usePrevSol Evaluate function with solution of current or previous time step
-     */
-    void computeStorage(PrimaryVariables &storage, const int scvIdx, bool usePrevSol) const
-    {
-        // compute the storage term for phase mass
-        ParentType::computeStorage(storage, scvIdx, usePrevSol);
-
-        // if flag usePrevSol is set, the solution from the previous
-        // time step is used, otherwise the current solution is
-        // used. The secondary variables are used accordingly.  This
-        // is required to compute the derivative of the storage term
-        // using the implicit euler method.
-        const ElementVolumeVariables &elemVolVars = usePrevSol ? this->prevVolVars_() : this->curVolVars_();
-        const VolumeVariables &volVars = elemVolVars[scvIdx];
-
-        // compute the energy storage
-        storage[energyEqIdx] = volVars.porosity()
-            *(
-                volVars.density(wPhaseIdx)
-                *volVars.internalEnergy(wPhaseIdx)
-                *volVars.saturation(wPhaseIdx)
-                +
-                volVars.density(nPhaseIdx)
-                *volVars.internalEnergy(nPhaseIdx)
-                *volVars.saturation(nPhaseIdx)
-                +
-                volVars.density(gPhaseIdx)
-                *volVars.internalEnergy(gPhaseIdx)
-                *volVars.saturation(gPhaseIdx)
-            )
-            + volVars.temperature()*volVars.heatCapacity();
-    }
-
-    /*!
-     * \brief Evaluates the advective mass flux and the heat flux
-     *        over a face of a subcontrol volume and writes the result in
-     *        the flux vector.
-     *
-     * \param flux The advective flux over the SCV (sub-control-volume) face for each component
-     * \param fluxVars The flux variables at the current SCV face
-     *
-     * This method is called by compute flux (base class).
-     */
-    void computeAdvectiveFlux(PrimaryVariables &flux,
-                              const FluxVariables &fluxVars) const
-    {
-        // advective mass flux
-        ParentType::computeAdvectiveFlux(flux, fluxVars);
-
-        static const Scalar massUpwindWeight = GET_PARAM_FROM_GROUP(TypeTag, Scalar, Implicit, MassUpwindWeight);
-
-        // advective heat flux in all phases
-        flux[energyEqIdx] = 0;
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-            // vertex data of the upstream and the downstream vertices
-            const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(phaseIdx));
-            const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(phaseIdx));
-
-            flux[energyEqIdx] +=
-                fluxVars.volumeFlux(phaseIdx) * (
-                                              massUpwindWeight * // upstream vertex
-                                              (  up.density(phaseIdx) *
-                                                 up.enthalpy(phaseIdx))
-                                              +
-                                              (1-massUpwindWeight) * // downstream vertex
-                                              (  dn.density(phaseIdx) *
-                                                 dn.enthalpy(phaseIdx)) );
-        }
-    }
-
-    /*!
-     * \brief Adds the diffusive heat flux to the flux vector over
-     *        the face of a sub-control volume.
-     *
-     * \param flux The diffusive flux over the sub-control-volume face for each conservation quantity (mass, energy)
-     * \param fluxVars The flux variables at the current sub-control-volume face
-     *
-     * This method is called by compute flux (base class).
-     */
-    void computeDiffusiveFlux(PrimaryVariables &flux,
-                              const FluxVariables &fluxVars) const
-    {
-        // diffusive mass flux
-        ParentType::computeDiffusiveFlux(flux, fluxVars);
-
-        // diffusive heat flux
-        flux[temperatureIdx] +=
-            fluxVars.normalMatrixHeatFlux();
-    }
-};
-
-}
-
-#endif
diff --git a/dumux/implicit/3p3cni/3p3cnimodel.hh b/dumux/implicit/3p3cni/3p3cnimodel.hh
deleted file mode 100644
index 9bf5f95d70e3e30588034ca2fb1f7031683cdb15..0000000000000000000000000000000000000000
--- a/dumux/implicit/3p3cni/3p3cnimodel.hh
+++ /dev/null
@@ -1,112 +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 Adaption of the fully implicit scheme to the non-isothermal
- *        three-phase three-component flow model.
- */
-#ifndef DUMUX_NEW_3P3CNI_MODEL_HH
-#define DUMUX_NEW_3P3CNI_MODEL_HH
-
-#include "3p3cniproperties.hh"
-#include <dumux/implicit/3p3c/3p3cmodel.hh>
-
-namespace Dumux {
-/*!
- * \ingroup ThreePThreeCNIModel
- * \brief Adaption of the fully implicit scheme to the non-isothermal
- *        three-phase three-component flow model.
- *
- * This model implements three-phase three-component flow of three fluid phases
- * \f$\alpha \in \{ \text{water, gas, NAPL} \}\f$ each composed of up to three components
- * \f$\kappa \in \{ \text{water, air, contaminant} \}\f$. The standard multiphase Darcy
- * approach is used as the equation for the conservation of momentum:
- * \f[
- v_\alpha = - \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K}
- \left(\textbf{grad}\, p_\alpha - \varrho_{\alpha} \mbox{\bf g} \right)
- * \f]
- *
- * By inserting this into the equations for the conservation of the
- * components, one transport equation for each component is obtained as
- * \f{eqnarray*}
- && \phi \frac{\partial (\sum_\alpha \varrho_\alpha X_\alpha^\kappa
- S_\alpha )}{\partial t}
- - \sum\limits_\alpha \text{div} \left\{ \frac{k_{r\alpha}}{\mu_\alpha}
- \varrho_\alpha X_\alpha^\kappa \mathbf{K}
- (\textbf{grad}\; p_\alpha - \varrho_\alpha \mbox{\bf g}) \right\}
- \nonumber \\
- \nonumber \\
- && - \sum\limits_\alpha \text{div} \left\{ D_{\alpha,\text{pm}}^\kappa
- \varrho_\alpha \frac{M^\kappa}{M_\alpha} \textbf{grad} x^\kappa_{\alpha} \right\}
- - q^\kappa = 0 \qquad \forall \kappa , \; \forall \alpha
- \f}
- *
- * Note that these balance equations above are molar.
- * In addition to that, a single balance of thermal energy is formulated
- * for the fluid-filled porous medium under the assumption of local thermal
- * equilibrium
- * \f{eqnarray*}
- && \phi \frac{\partial \left( \sum_\alpha \varrho_\alpha u_\alpha S_\alpha \right)}{\partial t}
- + \left( 1 - \phi \right) \frac{\partial (\varrho_s c_s T)}{\partial t}
- - \sum_\alpha \text{div} \left\{ \varrho_\alpha h_\alpha
- \frac{k_{r\alpha}}{\mu_\alpha} \mathbf{K} \left( \textbf{grad}\,
- p_\alpha
- - \varrho_\alpha \mathbf{g} \right) \right\} \\
- &-& \text{div} \left( \lambda_{pm} \textbf{grad} \, T \right)
- - q^h = 0 \qquad \alpha \in \{w, n, g\}
- \f}
- *
- * All equations are discretized using a vertex-centered finite volume (box)
- * or cell-centered finite volume scheme as spatial
- * and the implicit Euler method as time discretization.
- *
- * The model uses commonly applied auxiliary conditions like
- * \f$S_w + S_n + S_g = 1\f$ for the saturations and
- * \f$x^w_\alpha + x^a_\alpha + x^c_\alpha = 1\f$ for the mole fractions.
- * Furthermore, the phase pressures are related to each other via
- * capillary pressures between the fluid phases, which are functions of
- * the saturation, e.g. according to the approach of Parker et al.
- *
- * The used primary variables are dependent on the locally present fluid phases.
- * An adaptive primary variable switch is included. The phase state is stored for all nodes
- * of the system. The following cases can be distinguished:
- * <ul>
- *  <li> All three phases are present: Primary variables are two saturations \f$(S_w\f$ and \f$S_n)\f$,
- *       a pressure, in this case \f$p_g\f$, and the temperature \f$T\f$. </li>
- *  <li> Only the water phase is present: Primary variables are now the mole fractions of air and
- *       contaminant in the water phase \f$(x_w^a\f$ and \f$x_w^c)\f$, as well as temperature and the gas pressure,
- *       which is, of course, in a case where only the water phase is present, just the same as the water pressure. </li>
- *  <li> Gas and NAPL phases are present: Primary variables \f$(S_n\f$, \f$x_g^w\f$, \f$p_g\f$, \f$T)\f$. </li>
- *  <li> Water and NAPL phases are present: Primary variables \f$(S_n\f$, \f$x_w^a\f$, \f$p_g\f$, \f$T)\f$. </li>
- *  <li> Only gas phase is present: Primary variables \f$(x_g^w\f$, \f$x_g^c\f$, \f$p_g\f$, \f$T)\f$. </li>
- *  <li> Water and gas phases are present: Primary variables \f$(S_w\f$, \f$x_w^g\f$, \f$p_g\f$, \f$T)\f$. </li>
- * </ul>
- *
- */
-template<class TypeTag>
-class ThreePThreeCNIModel : public ThreePThreeCModel<TypeTag>
-{
-};
-
-}
-
-#include "3p3cnipropertydefaults.hh"
-
-#endif
diff --git a/dumux/implicit/3p3cni/3p3cniproperties.hh b/dumux/implicit/3p3cni/3p3cniproperties.hh
deleted file mode 100644
index 44ddd5df7e840799086c928af325e147b6e12685..0000000000000000000000000000000000000000
--- a/dumux/implicit/3p3cni/3p3cniproperties.hh
+++ /dev/null
@@ -1,51 +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/>.   *
- *****************************************************************************/
-/*!
- * \ingroup Properties
- * \ingroup ImplicitProperties
- * \ingroup ThreePThreeCNIModel
- * \file
- *
- * \brief Defines the properties required for the non-isothermal three-phase
- *        three-component fully implicit model.
- */
-#ifndef DUMUX_3P3CNI_PROPERTIES_HH
-#define DUMUX_3P3CNI_PROPERTIES_HH
-
-#warning You are including the old 2pni model which will be removed after 2.6. See CHANGELOG for details.
-
-#include <dumux/implicit/3p3c/3p3cproperties.hh>
-
-namespace Dumux
-{
-
-namespace Properties
-{
-//////////////////////////////////////////////////////////////////
-// Type tags
-//////////////////////////////////////////////////////////////////
-
-//! The type tags for the implicit non-isothermal three-phase three-component problems
-NEW_TYPE_TAG(ThreePThreeCNI, INHERITS_FROM(ThreePThreeC));
-NEW_TYPE_TAG(BoxThreePThreeCNI, INHERITS_FROM(BoxModel, ThreePThreeCNI));
-NEW_TYPE_TAG(CCThreePThreeCNI, INHERITS_FROM(CCModel, ThreePThreeCNI));
-}
-}
-
-#endif
diff --git a/dumux/implicit/3p3cni/3p3cnipropertydefaults.hh b/dumux/implicit/3p3cni/3p3cnipropertydefaults.hh
deleted file mode 100644
index 7fc8d9cb824d02b80f07e67217017447b8103329..0000000000000000000000000000000000000000
--- a/dumux/implicit/3p3cni/3p3cnipropertydefaults.hh
+++ /dev/null
@@ -1,73 +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/>.   *
- *****************************************************************************/
-/*!
- * \ingroup Properties
- * \ingroup ImplicitProperties
- * \ingroup ThreePThreeCNIModel
- */
-/*!
- * \file
- *
- * \brief Defines default values for most properties required by the
- *        non-isothermal three-phase three-component fully implicit model.
- */
-#ifndef DUMUX_3P3CNI_PROPERTY_DEFAULTS_HH
-#define DUMUX_3P3CNI_PROPERTY_DEFAULTS_HH
-
-#include <dumux/implicit/3p/3ppropertydefaults.hh>
-
-#include "3p3cnimodel.hh"
-#include "3p3cniindices.hh"
-#include "3p3cnilocalresidual.hh"
-#include "3p3cnivolumevariables.hh"
-#include "3p3cnifluxvariables.hh"
-
-namespace Dumux
-{
-
-namespace Properties
-{
-//////////////////////////////////////////////////////////////////
-// Property values
-//////////////////////////////////////////////////////////////////
-
-SET_INT_PROP(ThreePThreeCNI, NumEq, 4); //!< set the number of equations to 4
-
-//! Use the 3p3cni local jacobian operator for the 3p3cni model
-SET_TYPE_PROP(ThreePThreeCNI,
-              LocalResidual,
-              ThreePThreeCNILocalResidual<TypeTag>);
-
-//! the Model property
-SET_TYPE_PROP(ThreePThreeCNI, Model, ThreePThreeCNIModel<TypeTag>);
-
-//! the VolumeVariables property
-SET_TYPE_PROP(ThreePThreeCNI, VolumeVariables, ThreePThreeCNIVolumeVariables<TypeTag>);
-
-
-//! the FluxVariables property
-SET_TYPE_PROP(ThreePThreeCNI, FluxVariables, ThreePThreeCNIFluxVariables<TypeTag>);
-
-//! The indices required by the non-isothermal 3p3c model
-SET_TYPE_PROP(ThreePThreeCNI, Indices, ThreePThreeCNIIndices<TypeTag, 0>);
-
-}
-
-}
-#endif
diff --git a/dumux/implicit/3p3cni/3p3cnivolumevariables.hh b/dumux/implicit/3p3cni/3p3cnivolumevariables.hh
deleted file mode 100644
index 17a79cdfb3a598e124be3e402c59cd65adb3baa4..0000000000000000000000000000000000000000
--- a/dumux/implicit/3p3cni/3p3cnivolumevariables.hh
+++ /dev/null
@@ -1,158 +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 Contains the quantities which are constant within a
- *        finite volume in the non-isothermal three-phase three-component
- *        model.
- */
-#ifndef DUMUX_3P3CNI_VOLUME_VARIABLES_HH
-#define DUMUX_3P3CNI_VOLUME_VARIABLES_HH
-
-#include <dumux/implicit/3p3c/3p3cvolumevariables.hh>
-
-
-namespace Dumux
-{
-
-/*!
- * \ingroup ThreePThreeCNIModel
- * \ingroup ImplicitVolumeVariables
- * \brief Contains the quantities which are are constant within a
- *        finite volume in the non-isothermal three-phase three-component
- *        model.
- */
-template <class TypeTag>
-class ThreePThreeCNIVolumeVariables : public ThreePThreeCVolumeVariables<TypeTag>
-{
-    //! \cond false
-    typedef ThreePThreeCVolumeVariables<TypeTag> ParentType;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-
-    typedef typename GridView::template Codim<0>::Entity Element;
-    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-
-    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
-    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases) };
-    enum { temperatureIdx = Indices::temperatureIdx };
-
-    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
-    //! \endcond
-
-public:
-    /*!
-     * \copydoc ImplicitVolumeVariables::update
-     */
-    void update(const PrimaryVariables &priVars,
-                const Problem &problem,
-                const Element &element,
-                const FVElementGeometry &fvGeometry,
-                const int scvIdx,
-                bool isOldSol)
-    {
-        // vertex update data for the mass balance
-        ParentType::update(priVars,
-                           problem,
-                           element,
-                           fvGeometry,
-                           scvIdx,
-                           isOldSol);
-
-        typename FluidSystem::ParameterCache paramCache;
-        paramCache.updateAll(this->fluidState());
-
-        // the internal energies and the enthalpies
-        for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-            Scalar h = FluidSystem::enthalpy(this->fluidState_, paramCache, phaseIdx);
-            this->fluidState_.setEnthalpy(phaseIdx, h);
-        }
-    };
-
-    /*!
-     * \brief Returns the total internal energy of a phase in the
-     *        sub-control volume.
-     *
-     * \param phaseIdx The phase index
-     */
-    Scalar internalEnergy(int phaseIdx) const
-    { return this->fluidState_.internalEnergy(phaseIdx); };
-
-    /*!
-     * \brief Returns the total enthalpy of a phase in the sub-control
-     *        volume.
-     *
-     * \param phaseIdx The phase index
-     */
-    Scalar enthalpy(int phaseIdx) const
-    { return this->fluidState_.enthalpy(phaseIdx); };
-
-    /*!
-     * \brief Returns the total heat capacity \f$\mathrm{[J/(K*m^3]}\f$ of the rock matrix in
-     *        the sub-control volume.
-     */
-    Scalar heatCapacity() const
-    { return heatCapacity_; };
-
-protected:
-    // this method gets called by the parent class. since this method
-    // is protected, we are friends with our parent..
-    friend class ThreePThreeCVolumeVariables<TypeTag>;
-
-    static Scalar temperature_(const PrimaryVariables &priVars,
-                               const Problem& problem,
-                               const Element &element,
-                               const FVElementGeometry &fvGeometry,
-                               const int scvIdx)
-    {
-        return priVars[temperatureIdx];
-    }
-
-    /*!
-     * \brief Update all quantities for a given control volume.
-     *
-     * \param priVars The solution primary variables
-     * \param problem The problem
-     * \param element The element
-     * \param fvGeometry Evaluate function with solution of current or previous time step
-     * \param scvIdx The local index of the SCV (sub-control volume)
-     * \param isOldSol Evaluate function with solution of current or previous time step
-     */
-    void updateEnergy_(const PrimaryVariables &priVars,
-                       const Problem &problem,
-                       const Element &element,
-                       const FVElementGeometry &fvGeometry,
-                       const int scvIdx,
-                       bool isOldSol)
-    {
-        // copmute and set the heat capacity of the solid phase
-        heatCapacity_ = problem.spatialParams().heatCapacity(element, fvGeometry, scvIdx);
-        Valgrind::CheckDefined(heatCapacity_);
-    };
-
-    Scalar heatCapacity_;
-};
-
-} // end namespace
-
-#endif
diff --git a/dumux/implicit/3p3cni/Makefile.am b/dumux/implicit/3p3cni/Makefile.am
deleted file mode 100644
index c960520f8b1db14ab45fec83bac3e7fca8cf5c81..0000000000000000000000000000000000000000
--- a/dumux/implicit/3p3cni/Makefile.am
+++ /dev/null
@@ -1,4 +0,0 @@
-3p3cnidir = $(includedir)/dumux/implicit/3p3cni
-3p3cni_HEADERS := $(wildcard *.hh)
-
-include $(top_srcdir)/am/global-rules
diff --git a/dumux/implicit/co2ni/Makefile.am b/dumux/implicit/co2ni/Makefile.am
deleted file mode 100644
index 84015b26f806823ba88bba0e1d8283110a9c55a9..0000000000000000000000000000000000000000
--- a/dumux/implicit/co2ni/Makefile.am
+++ /dev/null
@@ -1,4 +0,0 @@
-co2nidir = $(includedir)/dumux/implicit/co2ni
-co2ni_HEADERS := $(wildcard *.hh)
-
-include $(top_srcdir)/am/global-rules
diff --git a/dumux/implicit/co2ni/co2nimodel.hh b/dumux/implicit/co2ni/co2nimodel.hh
deleted file mode 100644
index 3c32d9427b36900d8f8888d9e5dd32c5f6780b90..0000000000000000000000000000000000000000
--- a/dumux/implicit/co2ni/co2nimodel.hh
+++ /dev/null
@@ -1,59 +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 Adaption of the fully implicit scheme to the non-isothermal
- *        CO2 model.
- */
-#ifndef DUMUX_CO2NI_MODEL_HH
-#define DUMUX_CO2NI_MODEL_HH
-
-#include <dumux/implicit/2p2cni/2p2cniproperties.hh>
-#include <dumux/implicit/co2/co2model.hh>
-
-#warning You are including the old CO2NI model which will be removed after 2.6. See CHANGELOG for details.
-
-namespace Dumux {
-/*!
- * \ingroup CO2NIModel
-* \brief Adaption of the fully implicit scheme to the non-isothermal
- *        CO2 model.
- *
- *   See TwoPTwoCNI model for reference to the equations.
- *   The CO2NI model is derived from the CO2 model. In the CO2 model the phase switch criterion
- *   is different from the 2p2c model.
- *   The phase switch occurs when the equilibrium concentration
- *   of a component in a phase is exceeded instead of the sum of the components in the virtual phase
- *   (the phase which is not present) being greater that unity as done in the 2p2c model.
- *   The CO2VolumeVariables do not use a constraint solver for calculating the mole fractions as is the
- *   case in the 2p2c model. Instead mole fractions are calculated in the FluidSystem with a given
- *   temperature, pressure and salinity.
- *
- */
-template<class TypeTag>
-class CO2NIModel : public CO2Model<TypeTag>
-{
-};
-
-}
-
-#include <dumux/implicit/2p2cni/2p2cnipropertydefaults.hh>
-
-#endif
diff --git a/dumux/implicit/co2ni/co2nivolumevariables.hh b/dumux/implicit/co2ni/co2nivolumevariables.hh
deleted file mode 100644
index 35f8420b131c82e58d111cb32d28c199471ed938..0000000000000000000000000000000000000000
--- a/dumux/implicit/co2ni/co2nivolumevariables.hh
+++ /dev/null
@@ -1,138 +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 Contains the quantities which are constant within a
- *        finite volume in the non-isothermal CO2
- *        model.
- */
-#ifndef DUMUX_CO2NI_VOLUME_VARIABLES_HH
-#define DUMUX_CO2NI_VOLUME_VARIABLES_HH
-
-#include <dumux/implicit/2p2cni/2p2cniproperties.hh>
-#include <dumux/implicit/co2/co2volumevariables.hh>
-
-namespace Dumux
-{
-
-/*!
- * \ingroup CO2NIModel
- * \ingroup ImplicitVolumeVariables
- * \brief Contains the quantities which are are constant within a
- *        finite volume in the non-isothermal CO2
- *        model.
- */
-template <class TypeTag>
-class CO2NIVolumeVariables : public CO2VolumeVariables<TypeTag>
-{
-    //! \cond false
-    typedef CO2VolumeVariables<TypeTag> ParentType;
-    typedef typename ParentType::FluidState FluidState;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
-    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
-    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
-    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
-    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
-
-    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
-    enum { temperatureIdx = Indices::temperatureIdx };
-
-    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
-    typedef typename GridView::template Codim<0>::Entity Element;
-    //! \endcond
-
-public:
-    /*!
-     * \brief Returns the total internal energy of a phase in the
-     *        sub-control volume in \f$\mathrm{[J/kg]}\f$.
-     *
-     * \param phaseIdx The phase index
-     */
-    Scalar internalEnergy(const int phaseIdx) const
-    { return this->fluidState_.internalEnergy(phaseIdx); };
-
-    /*!
-     * \brief Returns the total enthalpy of a phase in the control
-     *        volume in \f$\mathrm{[J/kg]}\f$.
-     *
-     * \param phaseIdx The phase index
-     */
-    Scalar enthalpy(const int phaseIdx) const
-    { return this->fluidState_.enthalpy(phaseIdx); };
-
-    /*!
-     * \brief Returns the total heat capacity \f$\mathrm{[J/(K*m^3]}\f$ of the rock matrix in
-     *        the control volume.
-     */
-    Scalar heatCapacity() const
-    { return heatCapacity_; };
-
-    /*!
-     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m*K)]}\f$ of the fluid phase in
-     *        the control volume.
-     */
-    Scalar thermalConductivity(const int phaseIdx) const
-    { return FluidSystem::thermalConductivity(this->fluidState_, phaseIdx); };
-
-
-protected:
-    // this method gets called by the parent class. since this method
-    // is protected, we are friends with our parent..
-    friend class CO2VolumeVariables<TypeTag>;
-
-    static Scalar temperature_(const PrimaryVariables &priVars,
-                               const Problem& problem,
-                               const Element &element,
-                               const FVElementGeometry &fvGeometry,
-                               const int scvIdx)
-    {
-        return priVars[temperatureIdx];
-    }
-
-    template<class ParameterCache>
-    static Scalar enthalpy_(const FluidState& fluidState,
-                            const ParameterCache& paramCache,
-                            const int phaseIdx)
-    {
-        return FluidSystem::enthalpy(fluidState, paramCache, phaseIdx);
-    }
-
-    /*!
-     * \copydoc ImplicitVolumeVariables::update
-     */
-    void updateEnergy_(const PrimaryVariables &priVars,
-                       const Problem &problem,
-                       const Element &element,
-                       const FVElementGeometry &fvGeometry,
-                       const int scvIdx,
-                       bool isOldSol)
-    {
-        // copmute and set the heat capacity of the solid phase
-        heatCapacity_ = problem.spatialParams().heatCapacity(element, fvGeometry, scvIdx);
-        Valgrind::CheckDefined(heatCapacity_);
-    };
-
-    Scalar heatCapacity_;
-};
-
-} // end namespace
-
-#endif