diff --git a/dumux/porousmediumflow/richards/implicit/properties.hh b/dumux/porousmediumflow/richards/implicit/properties.hh
index 7c5184bc24d2dc821891515218c29da0c587f4bf..4ffaae0e2c67421c1abd10d785a4e7ce49b3adf8 100644
--- a/dumux/porousmediumflow/richards/implicit/properties.hh
+++ b/dumux/porousmediumflow/richards/implicit/properties.hh
@@ -148,13 +148,6 @@ public:
     using type = ImmiscibleFluidState<Scalar, FluidSystem>;
 };
 
-//! default value for the forchheimer coefficient
-// Source: Ward, J.C. 1964 Turbulent flow in porous media. ASCE J. Hydraul. Div 90.
-//        Actually the Forchheimer coefficient is also a function of the dimensions of the
-//        porous medium. Taking it as a constant is only a first approximation
-//        (Nield, Bejan, Convection in porous media, 2006, p. 10)
-SET_SCALAR_PROP(Richards, SpatialParamsForchCoeff, 0.55);
-
 //! Somerton is used as default model to compute the effective thermal heat conductivity
 SET_PROP(RichardsNI, ThermalConductivityModel)
 {
diff --git a/dumux/porousmediumflow/richardsnc/implicit/model.hh b/dumux/porousmediumflow/richardsnc/implicit/model.hh
index 51e576823cb1dbf94f665a36f58147446239a93c..2a0e950ffcd377e5ead3f2b5203059e045da0e3c 100644
--- a/dumux/porousmediumflow/richardsnc/implicit/model.hh
+++ b/dumux/porousmediumflow/richardsnc/implicit/model.hh
@@ -28,129 +28,4 @@
 
 #include "properties.hh"
 
-namespace Dumux
-{
-
-/*!
- * \ingroup RichardsModel
- *
- * \brief This model which implements a variant of the Richards'
- *        equation for quasi-twophase flow.
- *
- * In the unsaturated zone, Richards' equation
- \f[
- \frac{\partial\;\phi S_w \varrho_w}{\partial t}
- -
- \text{div} \left\lbrace
- \varrho_w \frac{k_{rw}}{\mu_w} \; \mathbf{K} \;
- \left( \text{\textbf{grad}}
- p_w - \varrho_w \textbf{g}
- \right)
- \right\rbrace
- =
- q_w,
- \f]
- * is frequently used to
- * approximate the water distribution above the groundwater level.
- *
- * It can be derived from the two-phase equations, i.e.
- \f[
- \phi\frac{\partial S_\alpha \varrho_\alpha}{\partial t}
- -
- \text{div} \left\lbrace
- \varrho_\alpha \frac{k_{r\alpha}}{\mu_\alpha}\; \mathbf{K} \;
- \left( \text{\textbf{grad}}
- p_\alpha - \varrho_\alpha \textbf{g}
- \right)
- \right\rbrace
- =
- q_\alpha,
- \f]
- * where \f$\alpha \in \{w, n\}\f$ is the fluid phase,
- * \f$\kappa \in \{ w, a \}\f$ are the components,
- * \f$\rho_\alpha\f$ is the fluid density, \f$S_\alpha\f$ is the fluid
- * saturation, \f$\phi\f$ is the porosity of the soil,
- * \f$k_{r\alpha}\f$ is the relative permeability for the fluid,
- * \f$\mu_\alpha\f$ is the fluid's dynamic viscosity, \f$\mathbf{K}\f$ is the
- * intrinsic permeability, \f$p_\alpha\f$ is the fluid pressure and
- * \f$g\f$ is the potential of the gravity field.
- *
- * In contrast to the full two-phase model, the Richards model assumes
- * gas as the non-wetting fluid and that it exhibits a much lower
- * viscosity than the (liquid) wetting phase. (For example at
- * atmospheric pressure and at room temperature, the viscosity of air
- * is only about \f$1\%\f$ of the viscosity of liquid water.) As a
- * consequence, the \f$\frac{k_{r\alpha}}{\mu_\alpha}\f$ term
- * typically is much larger for the gas phase than for the wetting
- * phase. For this reason, the Richards model assumes that
- * \f$\frac{k_{rn}}{\mu_n}\f$ is infinitly large. This implies that
- * the pressure of the gas phase is equivalent to the static pressure
- * distribution and that therefore, mass conservation only needs to be
- * considered for the wetting phase.
- *
- * The model thus choses the absolute pressure of the wetting phase
- * \f$p_w\f$ as its only primary variable. The wetting phase
- * saturation is calculated using the inverse of the capillary
- * pressure, i.e.
- \f[
- S_w = p_c^{-1}(p_n - p_w)
- \f]
- * holds, where \f$p_n\f$ is a given reference pressure. Nota bene,
- * that the last step is assumes that the capillary
- * pressure-saturation curve can be uniquely inverted, so it is not
- * possible to set the capillary pressure to zero when using the
- * Richards model!
- */
-
-template<class TypeTag >
-class RichardsNCModel : public GET_PROP_TYPE(TypeTag, BaseModel)
-{
-    using ParentType = typename GET_PROP_TYPE(TypeTag, BaseModel);
-    using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
-    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
-    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
-    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
-
-    static const int wPhaseIdx = Indices::wPhaseIdx;
-    static const int nPhaseIdx = Indices::nPhaseIdx;
-    static const int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
-
-public:
-
-    /*!
-     * \brief Apply the initial conditions to the model.
-     *
-     * \param problem The object representing the problem which needs to
-     *             be simulated.
-     */
-    void init(Problem& problem)
-    {
-        ParentType::init(problem);
-
-        // register standardized vtk output fields
-        auto& vtkOutputModule = problem.vtkOutputModule();
-        vtkOutputModule.addSecondaryVariable("Sw", [](const VolumeVariables& v){ return v.saturation(wPhaseIdx); });
-        vtkOutputModule.addSecondaryVariable("Sn", [](const VolumeVariables& v){ return v.saturation(nPhaseIdx); });
-        vtkOutputModule.addSecondaryVariable("pw", [](const VolumeVariables& v){ return v.pressure(wPhaseIdx); });
-        vtkOutputModule.addSecondaryVariable("pn", [](const VolumeVariables& v){ return v.pressure(nPhaseIdx); });
-        vtkOutputModule.addSecondaryVariable("pc", [](const VolumeVariables& v){ return v.capillaryPressure(); });
-        vtkOutputModule.addSecondaryVariable("density", [](const VolumeVariables& v){ return v.density(wPhaseIdx); });
-        vtkOutputModule.addSecondaryVariable("mobility", [](const VolumeVariables& v){ return v.mobility(wPhaseIdx); });
-        vtkOutputModule.addSecondaryVariable("kr", [](const VolumeVariables& v){ return v.relativePermeability(wPhaseIdx); });
-        vtkOutputModule.addSecondaryVariable("temperature", [](const VolumeVariables& v){ return v.temperature(); });
-        vtkOutputModule.addSecondaryVariable("porosity", [](const VolumeVariables& v){ return v.porosity(); });
-        if(GET_PARAM_FROM_GROUP(TypeTag, bool, Problem, EnableGravity));
-            vtkOutputModule.addSecondaryVariable("pressure head", [](const VolumeVariables& v){ return v.pressureHead(wPhaseIdx); });
-        vtkOutputModule.addSecondaryVariable("water content", [](const VolumeVariables& v){ return v.waterContent(wPhaseIdx); });
-
-        for (int k = 0; k < numComponents; ++k)
-            vtkOutputModule.addSecondaryVariable("x^" + FluidSystem::phaseName(wPhaseIdx) + "_" + FluidSystem::componentName(k),
-                                                 [k](const VolumeVariables& v){ return v.moleFraction(wPhaseIdx, k); });
-    }
-};
-
-} // end namespace Dumux
-
-#include "propertydefaults.hh"
-
 #endif
diff --git a/dumux/porousmediumflow/richardsnc/implicit/properties.hh b/dumux/porousmediumflow/richardsnc/implicit/properties.hh
index 4d3a5b76ff29c0390b927d4587182fc859d35879..eeb54907ef1457223ef72726baa49dd8f438b8b5 100644
--- a/dumux/porousmediumflow/richardsnc/implicit/properties.hh
+++ b/dumux/porousmediumflow/richardsnc/implicit/properties.hh
@@ -29,8 +29,27 @@
 #ifndef DUMUX_RICHARDSNC_PROPERTIES_HH
 #define DUMUX_RICHARDSNC_PROPERTIES_HH
 
+#include <dumux/common/basicproperties.hh>
+#include <dumux/linear/linearsolverproperties.hh>
+
+#include <dumux/porousmediumflow/compositional/localresidual.hh>
+#include <dumux/porousmediumflow/richards/implicit/newtoncontroller.hh>
+
+#include <dumux/material/spatialparams/implicit1p.hh>
+#include <dumux/material/fluidmatrixinteractions/diffusivitymillingtonquirk.hh>
+#include <dumux/material/fluidmatrixinteractions/1p/thermalconductivityaverage.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/components/constant.hh>
+#include <dumux/material/fluidsystems/liquidphase2c.hh>
+#include <dumux/material/fluidstates/compositional.hh>
+
+#include <dumux/porousmediumflow/properties.hh>
 #include <dumux/porousmediumflow/nonisothermal/implicit/properties.hh>
 
+#include "volumevariables.hh"
+#include "indices.hh"
+#include "vtkoutputfields.hh"
+
 namespace Dumux
 {
 // \{
@@ -42,34 +61,123 @@ namespace Properties
 //////////////////////////////////////////////////////////////////
 
 //! The type tags for the implicit isothermal one-phase two-component problems
-NEW_TYPE_TAG(RichardsNC);
-NEW_TYPE_TAG(BoxRichardsNC, INHERITS_FROM(BoxModel, RichardsNC));
-NEW_TYPE_TAG(CCRichardsNC, INHERITS_FROM(CCTpfaModel, RichardsNC));
+NEW_TYPE_TAG(RichardsNC, INHERITS_FROM(PorousMediumFlow, NumericModel, LinearSolverTypeTag));
+NEW_TYPE_TAG(RichardsNCNI, INHERITS_FROM(Richards, NonIsothermal));
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////
+// Property values
+//////////////////////////////////////////////////////////////////
+SET_PROP(RichardsNC, NumEq)
+{
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    static const int value = FluidSystem::numComponents;
+};
 
-//! The type tags for the corresponding non-isothermal problems
-NEW_TYPE_TAG(RichardsNCNI, INHERITS_FROM(RichardsNC, NonIsothermal));
-NEW_TYPE_TAG(BoxRichardsNCNI, INHERITS_FROM(BoxModel, RichardsNCNI));
-NEW_TYPE_TAG(CCRichardsNCNI, INHERITS_FROM(CCTpfaModel, RichardsNCNI));
+SET_PROP(RichardsNC, NumPhases)
+{
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    static const int value = FluidSystem::numPhases;
+    static_assert(value == 1, "You can only use one-phasic fluid systems with the Richards n-component model!");
+};
+
+SET_PROP(RichardsNC, NumComponents)
+{
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    static const int value = FluidSystem::numComponents;
+};
+
+SET_BOOL_PROP(RichardsNC, UseMoles, true); //!< Define that per default mole fractions are used in the balance equations
+
+//! Use the dedicated local residual
+SET_TYPE_PROP(RichardsNC, LocalResidual, CompositionalLocalResidual<TypeTag>);
+
+//! We set the replaceCompIdx to 0, i.e. the first equation is substituted with
+//! the total mass balance, i.e. the phase balance
+SET_INT_PROP(RichardsNC, ReplaceCompEqIdx, 0);
+
+//! define the VolumeVariables
+SET_TYPE_PROP(RichardsNC, VolumeVariables, RichardsNCVolumeVariables<TypeTag>);
+
+/*!
+ *\brief The fluid system used by the model.
+ *
+ * By default this uses the liquid phase fluid system with simple H2O.
+ */
+SET_PROP(RichardsNC, FluidSystem)
+{
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using type = FluidSystems::LiquidPhaseTwoC<TypeTag, SimpleH2O<Scalar>, Constant<TypeTag, Scalar>>;
+};
+
+/*!
+ * \brief The fluid state which is used by the volume variables to
+ *        store the thermodynamic state. This should be chosen
+ *        appropriately for the model ((non-)isothermal, equilibrium, ...).
+ *        This can be done in the problem.
+ */
+SET_PROP(RichardsNC, FluidState)
+{
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    using type = CompositionalFluidState<Scalar, FluidSystem>;
+};
+SET_TYPE_PROP(RichardsNC, VtkOutputFields, RichardsNCVtkOutputFields<TypeTag>);           //! Set the vtk output fields specific to the twop model
+/*!
+ * \brief Set type of the parameter objects for the material law
+ *
+ * By default this is just retrieved from the material law.
+ */
+SET_TYPE_PROP(RichardsNC, MaterialLawParams, typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params);
+
+//! Set the indices used
+SET_TYPE_PROP(RichardsNC, Indices, RichardsNCIndices<TypeTag>);
+//! The spatial parameters to be employed.
+//! Use ImplicitSpatialParamsOneP by default.
+SET_TYPE_PROP(RichardsNC, SpatialParams, ImplicitSpatialParamsOneP<TypeTag>);
+
+//! The model after Millington (1961) is used for the effective diffusivity
+SET_PROP(RichardsNC, EffectiveDiffusivityModel)
+{
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using type = DiffusivityMillingtonQuirk<Scalar>;
+};
+
+//physical processes to be considered by the isothermal model
+SET_BOOL_PROP(RichardsNC, EnableAdvection, true);
+SET_BOOL_PROP(RichardsNC, EnableMolecularDiffusion, true);
+SET_BOOL_PROP(RichardsNC, EnableEnergyBalance, false);
+
+//! average is used as default model to compute the effective thermal heat conductivity
+SET_PROP(RichardsNCNI, ThermalConductivityModel)
+{
+    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+    using type = ThermalConductivityAverage<Scalar>;
+};
 
 //////////////////////////////////////////////////////////////////
-// Property tags
+// Property values for isothermal model required for the general non-isothermal model
 //////////////////////////////////////////////////////////////////
 
-NEW_PROP_TAG(NumPhases);   //!< Number of fluid phases in the system
-NEW_PROP_TAG(NumComponents);   //!< Number of fluid components in the system
-NEW_PROP_TAG(MaterialLaw);   //!< The material law which ought to be used (by default extracted from the spatial parameters)
-NEW_PROP_TAG(Indices); //!< Enumerations for the model
-NEW_PROP_TAG(SpatialParams); //!< The type of the spatial parameters
-NEW_PROP_TAG(EffectiveDiffusivityModel); //!< The employed model for the computation of the effective diffusivity
-NEW_PROP_TAG(FluidSystem); //!< Type of the multi-component relations
-NEW_PROP_TAG(FluidState); //!< Type of the fluid state to be used
-NEW_PROP_TAG(ProblemEnableGravity); //!< Returns whether gravity is considered in the problem
-NEW_PROP_TAG(UseMoles); //!< Defines whether mole (true) or mass (false) fractions are used
-NEW_PROP_TAG(SpatialParamsForchCoeff); //!< Property for the forchheimer coefficient
-NEW_PROP_TAG(TauTortuosity); //!< Tortuosity value (tau) used in macroscopic diffusion
-
-}
+//set isothermal VolumeVariables
+SET_TYPE_PROP(RichardsNCNI, IsothermalVolumeVariables, RichardsNCVolumeVariables<TypeTag>);
+
+//set isothermal LocalResidual
+SET_TYPE_PROP(RichardsNCNI, IsothermalLocalResidual, CompositionalLocalResidual<TypeTag>);
+
+//set isothermal Indices
+SET_TYPE_PROP(RichardsNCNI, IsothermalIndices, RichardsNCIndices<TypeTag>);
+
+//set isothermal NumEq
+SET_PROP(RichardsNCNI, IsothermalNumEq)
+{
+    using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+    static const int value = FluidSystem::numComponents;
+};
+
+} // end namespace Properties
 // \}
-}
+} // end namespace Dumux
 
 #endif