diff --git a/dumux/porousmediumflow/CMakeLists.txt b/dumux/porousmediumflow/CMakeLists.txt
index 76eeae0efb74a908bf13f572a44bea7a332b2673..ff813a915a6f9c2750cc6438dd91feb96c72473d 100644
--- a/dumux/porousmediumflow/CMakeLists.txt
+++ b/dumux/porousmediumflow/CMakeLists.txt
@@ -16,5 +16,6 @@ add_subdirectory("implicit")
add_subdirectory("mpnc")
add_subdirectory("nonisothermal")
add_subdirectory("richards")
+add_subdirectory("richardsnc")
add_subdirectory("sequential")
add_subdirectory("tracer")
diff --git a/dumux/porousmediumflow/richardsnc/CMakeLists.txt b/dumux/porousmediumflow/richardsnc/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ba8341c614f1a2c797c95f5402f602025f1087b1
--- /dev/null
+++ b/dumux/porousmediumflow/richardsnc/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory("implicit")
diff --git a/dumux/porousmediumflow/richardsnc/implicit/CMakeLists.txt b/dumux/porousmediumflow/richardsnc/implicit/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c19a6e9a869fda004b50d03b5387123ad6f422f4
--- /dev/null
+++ b/dumux/porousmediumflow/richardsnc/implicit/CMakeLists.txt
@@ -0,0 +1,10 @@
+
+#install headers
+install(FILES
+indices.hh
+model.hh
+localresidual.hh
+properties.hh
+propertydefaults.hh
+volumevariables.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/porousmediumflow/richardsnc/implicit)
diff --git a/dumux/porousmediumflow/richardsnc/implicit/indices.hh b/dumux/porousmediumflow/richardsnc/implicit/indices.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f3912564cc909f7e508aa7e0b2ee43d43014ffea
--- /dev/null
+++ b/dumux/porousmediumflow/richardsnc/implicit/indices.hh
@@ -0,0 +1,63 @@
+// -*- 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 primary variable and equation indices used by
+ * the richardsnc model
+ */
+
+#ifndef DUMUX_RICHARDSNC_INDICES_HH
+#define DUMUX_RICHARDSNC_INDICES_HH
+
+#include "properties.hh"
+
+namespace Dumux
+{
+// \{
+
+/*!
+ * \ingroup RichardsNCModel
+ * \ingroup ImplicitIndices
+ * \brief The indices for the isothermal Richards, n-component model.
+ */
+template <class TypeTag, int PVOffset = 0>
+struct RichardsNCIndices
+{
+
+ //! Set the index of the phases for accessing the volvars
+ static const int wPhaseIdx = 0;
+ static const int nPhaseIdx = 1;
+
+ //! Component indices
+ static const int compMainIdx = PVOffset + 0; //!< main component index
+
+ //! primary variable indices
+ static const int pressureIdx = PVOffset + 0; //!< pressure
+
+ //! \note These indices make sense if the first balance is replaced by the
+ //! total mass balance.
+
+ //! Equation indices
+ static const int conti0EqIdx = PVOffset + 0; //!< continuity equation index
+};
+
+// \}
+}
+
+#endif
diff --git a/dumux/porousmediumflow/richardsnc/implicit/model.hh b/dumux/porousmediumflow/richardsnc/implicit/model.hh
new file mode 100644
index 0000000000000000000000000000000000000000..51e576823cb1dbf94f665a36f58147446239a93c
--- /dev/null
+++ b/dumux/porousmediumflow/richardsnc/implicit/model.hh
@@ -0,0 +1,156 @@
+// -*- 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 Base class for all models which use the Richards,
+ * n-component fully implicit model.
+ */
+
+#ifndef DUMUX_RICHARDSNC_MODEL_HH
+#define DUMUX_RICHARDSNC_MODEL_HH
+
+#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
new file mode 100644
index 0000000000000000000000000000000000000000..386bb184a65203a80cde087bacd68d6b514cfd43
--- /dev/null
+++ b/dumux/porousmediumflow/richardsnc/implicit/properties.hh
@@ -0,0 +1,79 @@
+// -*- 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 RichardsNCModel
+ * \file
+ *
+ * \brief Defines the properties required for the Richards,
+ * n-component fully implicit model.
+ */
+
+#ifndef DUMUX_RICHARDSNC_PROPERTIES_HH
+#define DUMUX_RICHARDSNC_PROPERTIES_HH
+
+
+#include <dumux/implicit/box/properties.hh>
+#include <dumux/implicit/cellcentered/tpfa/properties.hh>
+#include <dumux/porousmediumflow/nonisothermal/implicit/properties.hh>
+
+namespace Dumux
+{
+// \{
+namespace Properties
+{
+
+//////////////////////////////////////////////////////////////////
+// Type tags
+//////////////////////////////////////////////////////////////////
+
+//! 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));
+
+//! 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));
+
+//////////////////////////////////////////////////////////////////
+// Property tags
+//////////////////////////////////////////////////////////////////
+
+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(MaterialLawParams); //!< The type of the parameter object for the material law (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
+
+}
+// \}
+}
+
+#endif
diff --git a/dumux/porousmediumflow/richardsnc/implicit/propertydefaults.hh b/dumux/porousmediumflow/richardsnc/implicit/propertydefaults.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e5f51c13615a09382de584b5c50db39c15062c7b
--- /dev/null
+++ b/dumux/porousmediumflow/richardsnc/implicit/propertydefaults.hh
@@ -0,0 +1,194 @@
+// -*- 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 RichardsNCModel
+ * \file
+ *
+ * \brief Defines some default values for the properties of the
+ * Richards, n-component fully implicit model.
+ */
+
+#ifndef DUMUX_RICHARDSNC_PROPERTY_DEFAULTS_HH
+#define DUMUX_RICHARDSNC_PROPERTY_DEFAULTS_HH
+
+#include "properties.hh"
+#include "model.hh"
+#include "volumevariables.hh"
+#include "indices.hh"
+
+#include <dumux/porousmediumflow/compositional/localresidual.hh>
+#include <dumux/porousmediumflow/richards/implicit/newtoncontroller.hh>
+#include <dumux/porousmediumflow/nonisothermal/implicit/propertydefaults.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>
+
+namespace Dumux
+{
+// \{
+namespace Properties
+{
+//////////////////////////////////////////////////////////////////
+// Property values
+//////////////////////////////////////////////////////////////////
+SET_PROP(RichardsNC, NumEq)
+{
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ static const int value = FluidSystem::numComponents;
+};
+
+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 model
+SET_TYPE_PROP(RichardsNC, Model, RichardsNCModel<TypeTag>);
+
+//! define the VolumeVariables
+SET_TYPE_PROP(RichardsNC, VolumeVariables, RichardsNCVolumeVariables<TypeTag>);
+
+//! Smarter newton controller
+SET_TYPE_PROP(RichardsNC, NewtonController, RichardsNewtonController<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>;
+};
+
+/*!
+ * \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>;
+};
+
+// enable gravity by default
+SET_BOOL_PROP(RichardsNC, ProblemEnableGravity, true);
+
+//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);
+
+//! 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(RichardsNC, SpatialParamsForchCoeff, 0.55);
+
+/*!
+ * \brief default value for tortuosity value (tau) used in macroscopic diffusion
+ *
+ * Value is 0.5 according to Carman 1937: <i>Fluid flow through granular beds</i>
+ * \cite carman1937
+ */
+SET_SCALAR_PROP(RichardsNC, TauTortuosity, 0.5);
+
+//! 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 values for isothermal model required for the general non-isothermal model
+//////////////////////////////////////////////////////////////////
+
+// set isothermal Model
+SET_TYPE_PROP(RichardsNCNI, IsothermalModel, RichardsNCModel<TypeTag>);
+
+//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
diff --git a/dumux/porousmediumflow/richardsnc/implicit/volumevariables.hh b/dumux/porousmediumflow/richardsnc/implicit/volumevariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..e1f325215b7c9e1a4da9ca9d81c873de62e5c660
--- /dev/null
+++ b/dumux/porousmediumflow/richardsnc/implicit/volumevariables.hh
@@ -0,0 +1,392 @@
+// -*- 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 @copybrief Dumux::RichardsNCVolumeVariables
+ */
+#ifndef DUMUX_RICHARDSNC_VOLUME_VARIABLES_HH
+#define DUMUX_RICHARDSNC_VOLUME_VARIABLES_HH
+
+#include <dumux/discretization/volumevariables.hh>
+
+#include "properties.hh"
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup RichardsNCModel
+ * \ingroup ImplicitVolumeVariables
+ * \brief Contains the quantities which are constant within a
+ * finite volume in the Richards, n-component model.
+ */
+template <class TypeTag>
+class RichardsNCVolumeVariables : public ImplicitVolumeVariables<TypeTag>
+{
+ using ParentType = ImplicitVolumeVariables<TypeTag>;
+
+ using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
+ using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
+ using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
+ using Implementation = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
+ using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
+ using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
+ using ElementSolutionVector = typename GET_PROP_TYPE(TypeTag, ElementSolutionVector);
+ using SpatialParams = typename GET_PROP_TYPE(TypeTag, SpatialParams);
+ using PermeabilityType = typename SpatialParams::PermeabilityType;
+ using FluidSystem = typename GET_PROP_TYPE(TypeTag, FluidSystem);
+ using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
+ enum
+ {
+ wPhaseIdx = Indices::wPhaseIdx,
+ nPhaseIdx = Indices::nPhaseIdx,
+ pressureIdx = Indices::pressureIdx
+ };
+
+ static const bool useMoles = GET_PROP_VALUE(TypeTag, UseMoles);
+ static const int dim = GridView::dimension;
+ static const int dimWorld = GridView::dimensionworld;
+ static const int numComponents = GET_PROP_VALUE(TypeTag, NumComponents);
+
+ using MaterialLaw = typename GET_PROP_TYPE(TypeTag, MaterialLaw);
+ using GlobalPosition = Dune::FieldVector<Scalar,dimWorld>;
+ using Element = typename GridView::template Codim<0>::Entity;
+
+public:
+
+ using FluidState = typename GET_PROP_TYPE(TypeTag, FluidState);
+
+ /*!
+ * \copydoc ImplicitVolumeVariables::update
+ */
+ void update(const ElementSolutionVector &elemSol,
+ const Problem &problem,
+ const Element &element,
+ const SubControlVolume &scv)
+ {
+ ParentType::update(elemSol, problem, element, scv);
+
+ //calculate all secondary variables from the primary variables and store results in fluidstate
+ completeFluidState(elemSol, problem, element, scv, fluidState_);
+
+ const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
+ relativePermeabilityWetting_ = MaterialLaw::krw(materialParams, fluidState_.saturation(wPhaseIdx));
+
+ //! precompute the minimum capillary pressure (entry pressure)
+ minPc_ = MaterialLaw::pc(materialParams, /*Sw=*/ 1.0);
+ pn_ = problem.nonWettingReferencePressure();
+ porosity_ = problem.spatialParams().porosity(element, scv, elemSol);
+ permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
+
+ // Second instance of a parameter cache.
+ // Could be avoided if diffusion coefficients also
+ // became part of the fluid state.
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updatePhase(fluidState_, wPhaseIdx);
+
+ const int compIIdx = wPhaseIdx;
+ for (unsigned int compJIdx = 0; compJIdx < numComponents; ++compJIdx)
+ if(compIIdx != compJIdx)
+ setDiffusionCoefficient_(compJIdx,
+ FluidSystem::binaryDiffusionCoefficient(fluidState_,
+ paramCache,
+ wPhaseIdx,
+ compIIdx,
+ compJIdx));
+ }
+
+ /*!
+ * \copydoc ImplicitModel::completeFluidState
+ */
+ static void completeFluidState(const ElementSolutionVector &elemSol,
+ const Problem& problem,
+ const Element& element,
+ const SubControlVolume &scv,
+ FluidState& fluidState)
+ {
+ Scalar t = ParentType::temperature(elemSol, problem, element, scv);
+ fluidState.setTemperature(t);
+
+ const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
+ const auto& priVars = ParentType::extractDofPriVars(elemSol, scv);
+
+ // set the wetting pressure
+ fluidState.setPressure(wPhaseIdx, priVars[pressureIdx]);
+
+ // compute the capillary pressure to compute the saturation
+ // the material law takes care of cases where pc gets negative
+ const Scalar pc = problem.nonWettingReferencePressure() - fluidState.pressure(wPhaseIdx);
+ const Scalar sw = MaterialLaw::sw(materialParams, pc);
+ fluidState.setSaturation(wPhaseIdx, sw);
+
+ // set the mole/mass fractions
+ if(useMoles)
+ {
+ Scalar sumSecondaryFractions = 0.0;
+ for (int compIdx = 1; compIdx < numComponents; ++compIdx)
+ {
+ fluidState.setMoleFraction(wPhaseIdx, compIdx, priVars[compIdx]);
+ sumSecondaryFractions += priVars[compIdx];
+ }
+ fluidState.setMoleFraction(wPhaseIdx, 0, 1.0 - sumSecondaryFractions);
+ }
+ else
+ {
+ for (int compIdx = 1; compIdx < numComponents; ++compIdx)
+ fluidState.setMassFraction(wPhaseIdx, compIdx, priVars[compIdx]);
+ }
+
+ // density and viscosity
+ typename FluidSystem::ParameterCache paramCache;
+ paramCache.updateAll(fluidState);
+ fluidState.setDensity(wPhaseIdx, FluidSystem::density(fluidState, paramCache, wPhaseIdx));
+ fluidState.setViscosity(wPhaseIdx, FluidSystem::viscosity(fluidState, paramCache, wPhaseIdx));
+
+ // compute and set the enthalpy
+ fluidState.setEnthalpy(wPhaseIdx, Implementation::enthalpy(fluidState, paramCache, wPhaseIdx));
+ }
+
+ /*!
+ * \brief Return the fluid configuration at the given primary
+ * variables
+ */
+ const FluidState &fluidState() const
+ { return fluidState_; }
+
+ /*!
+ * \brief Return the temperature
+ */
+ Scalar temperature() const
+ { return fluidState_.temperature(); }
+
+ /*!
+ * \brief Returns the average porosity [] within the control volume.
+ *
+ * The porosity is defined as the ratio of the pore space to the
+ * total volume, i.e. \f[ \Phi := \frac{V_{pore}}{V_{pore} + V_{rock}} \f]
+ */
+ Scalar porosity() const
+ { return porosity_; }
+
+ /*!
+ * \brief Returns the permeability within the control volume in \f$[m^2]\f$.
+ */
+ const PermeabilityType& permeability() const
+ { return permeability_; }
+
+ /*!
+ * \brief Returns the average absolute saturation [] of a given
+ * fluid phase within the finite volume.
+ *
+ * The saturation of a fluid phase is defined as the fraction of
+ * the pore volume filled by it, i.e.
+ * \f[ S_\alpha := \frac{V_\alpha}{V_{pore}} = \phi \frac{V_\alpha}{V} \f]
+ *
+ * \param phaseIdx The index of the fluid phase
+ */
+ Scalar saturation(const int phaseIdx) const
+ { return phaseIdx == wPhaseIdx ? fluidState_.saturation(wPhaseIdx) : 1.0-fluidState_.saturation(wPhaseIdx); }
+
+ /*!
+ * \brief Returns the average mass density \f$\mathrm{[kg/m^3]}\f$ of a given
+ * fluid phase within the control volume.
+ *
+ * \param phaseIdx The index of the fluid phase
+ */
+ Scalar density(const int phaseIdx) const
+ { return phaseIdx == wPhaseIdx ? fluidState_.density(phaseIdx) : 0.0; }
+
+ /*!
+ * \brief Returns the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within
+ * the control volume.
+ *
+ * For the non-wetting phase (i.e. the gas phase), we assume
+ * infinite mobility, which implies that the non-wetting phase
+ * pressure is equal to the finite volume's reference pressure
+ * defined by the problem.
+ *
+ * \param phaseIdx The index of the fluid phase
+ */
+ Scalar pressure(const int phaseIdx) const
+ { return phaseIdx == wPhaseIdx ? fluidState_.pressure(phaseIdx) : pn_; }
+
+ /*!
+ * \brief Returns the effective mobility \f$\mathrm{[1/(Pa*s)]}\f$ of a given phase within
+ * the control volume.
+ *
+ * The mobility of a fluid phase is defined as the relative
+ * permeability of the phase (given by the chosen material law)
+ * divided by the dynamic viscosity of the fluid, i.e.
+ * \f[ \lambda_\alpha := \frac{k_{r\alpha}}{\mu_\alpha} \f]
+ *
+ * \param phaseIdx The index of the fluid phase
+ */
+ Scalar mobility(const int phaseIdx) const
+ { return relativePermeability(phaseIdx)/fluidState_.viscosity(phaseIdx); }
+
+ /*!
+ * \brief Returns the dynamic viscosity \f$\mathrm{[Pa*s]}\f$ of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The index of the fluid phase
+ * \note The non-wetting phase is infinitely mobile
+ */
+ Scalar viscosity(const int phaseIdx) const
+ { return phaseIdx == wPhaseIdx ? fluidState_.viscosity(wPhaseIdx) : 0.0; }
+
+ /*!
+ * \brief Returns relative permeability [-] of a given phase within
+ * the control volume.
+ *
+ * \param phaseIdx The index of the fluid phase
+ */
+ Scalar relativePermeability(const int phaseIdx) const
+ { return phaseIdx == wPhaseIdx ? relativePermeabilityWetting_ : 1.0; }
+
+ /*!
+ * \brief Returns the effective capillary pressure \f$\mathrm{[Pa]}\f$ within the
+ * control volume.
+ *
+ * The capillary pressure is defined as the difference in
+ * pressures of the non-wetting and the wetting phase, i.e.
+ * \f[ p_c = p_n - p_w \f]
+ *
+ * \note Capillary pressures are always larger than the entry pressure
+ * This regularization doesn't affect the residual in which pc is not needed.
+ */
+ Scalar capillaryPressure() const
+ {
+ using std::max;
+ return max(minPc_, pn_ - fluidState_.pressure(wPhaseIdx));
+ }
+
+ /*!
+ * \brief Returns the pressureHead \f$\mathrm{[cm]}\f$ of a given phase within
+ * the control volume.
+ *
+ * For the non-wetting phase (i.e. the gas phase), we assume
+ * infinite mobility, which implies that the non-wetting phase
+ * pressure is equal to the finite volume's reference pressure
+ * defined by the problem.
+ *
+ * \param phaseIdx The index of the fluid phase
+ * \note this function is here as a convenience to the user to not have to
+ * manually do a conversion. It is not correct if the density is not constant
+ * or the gravity different
+ */
+ Scalar pressureHead(const int phaseIdx) const
+ { return 100.0 *(pressure(phaseIdx) - pn_)/density(phaseIdx)/9.81; }
+
+ /*!
+ * \brief Returns the water content
+ * fluid phase within the finite volume.
+ *
+ * The water content is defined as the fraction of
+ * the saturation devided by the porosity
+
+ * \param phaseIdx The index of the fluid phase
+ * \note this function is here as a convenience to the user to not have to
+ * manually do a conversion.
+ */
+ Scalar waterContent(const int phaseIdx) const
+ { return saturation(phaseIdx) * porosity_; }
+
+ /*!
+ * \brief Return molar density \f$\mathrm{[mol/m^3]}\f$ the of the fluid phase.
+ *
+ * We always forward to the fluid state with the phaseIdx property (see class description).
+ */
+ Scalar molarDensity(const int phaseIdx) const
+ { return phaseIdx == wPhaseIdx ? fluidState_.molarDensity(phaseIdx) : 0.0; }
+
+ /*!
+ * \brief Return mole fraction \f$\mathrm{[mol/mol]}\f$ of a component in the phase.
+ *
+ * \param compIdx The index of the component
+ *
+ * We always forward to the fluid state with the phaseIdx property (see class description).
+ */
+ Scalar moleFraction(const int phaseIdx, const int compIdx) const
+ { return phaseIdx == wPhaseIdx ? fluidState_.moleFraction(phaseIdx, compIdx) : 0.0; }
+
+ /*!
+ * \brief Return mass fraction \f$\mathrm{[kg/kg]}\f$ of a component in the phase.
+ *
+ * \param compIdx The index of the component
+ *
+ * We always forward to the fluid state with the phaseIdx property (see class description).
+ */
+ Scalar massFraction(const int phaseIdx, const int compIdx) const
+ { return phaseIdx == wPhaseIdx ? fluidState_.massFraction(phaseIdx, compIdx) : 0.0; }
+
+ /*!
+ * \brief Return concentration \f$\mathrm{[mol/m^3]}\f$ of a component in the phase.
+ *
+ * \param compIdx The index of the component
+ *
+ * We always forward to the fluid state with the phaseIdx property (see class description).
+ */
+ Scalar molarity(const int phaseIdx, const int compIdx) const
+ { return phaseIdx == wPhaseIdx ? fluidState_.molarity(phaseIdx, compIdx) : 0.0; }
+
+ /*!
+ * \brief Return the binary diffusion coefficient \f$\mathrm{[m^2/s]}\f$ in the fluid.
+ */
+ Scalar diffusionCoefficient(const int phaseIdx, const int compIdx) const
+ {
+ assert(phaseIdx == wPhaseIdx);
+ assert(compIdx > wPhaseIdx);
+ return diffCoefficient_[compIdx-1];
+ }
+
+ /*!
+ * \brief Returns the dispersivity of the fluid's streamlines.
+ */
+ const GlobalPosition &dispersivity() const
+ { return dispersivity_; }
+
+protected:
+ FluidState fluidState_; //! the fluid state
+ Scalar relativePermeabilityWetting_; //! the relative permeability of the wetting phase
+ Scalar porosity_; //! the porosity
+ PermeabilityType permeability_; //! the instrinsic permeability
+ Scalar pn_; //! the reference non-wetting pressure
+ Scalar minPc_; //! the minimum capillary pressure (entry pressure)
+
+private:
+ void setDiffusionCoefficient_(int compIdx, Scalar d)
+ {
+ assert(compIdx > wPhaseIdx);
+ diffCoefficient_[compIdx-1] = d;
+ }
+
+ std::array<Scalar, numComponents-1> diffCoefficient_;
+ GlobalPosition dispersivity_;
+
+ Implementation &asImp_()
+ { return *static_cast<Implementation*>(this); }
+
+ const Implementation &asImp_() const
+ { return *static_cast<const Implementation*>(this); }
+};
+
+}// end namespace Dumux
+
+#endif
diff --git a/test/porousmediumflow/CMakeLists.txt b/test/porousmediumflow/CMakeLists.txt
index 8b6532dc0c4aac178b247a2f0566117cf352cc62..83c6b1991ced243d44a8272d142a0ed16d6f3feb 100644
--- a/test/porousmediumflow/CMakeLists.txt
+++ b/test/porousmediumflow/CMakeLists.txt
@@ -13,4 +13,5 @@ add_subdirectory("3pwateroil")
add_subdirectory("co2")
add_subdirectory("mpnc")
add_subdirectory("richards")
+add_subdirectory("richardsnc")
add_subdirectory("tracer")
diff --git a/test/porousmediumflow/richardsnc/CMakeLists.txt b/test/porousmediumflow/richardsnc/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ba8341c614f1a2c797c95f5402f602025f1087b1
--- /dev/null
+++ b/test/porousmediumflow/richardsnc/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory("implicit")
diff --git a/test/porousmediumflow/richardsnc/implicit/CMakeLists.txt b/test/porousmediumflow/richardsnc/implicit/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..8352fc02e6f0be792c0567b444b0f77d235ce16b
--- /dev/null
+++ b/test/porousmediumflow/richardsnc/implicit/CMakeLists.txt
@@ -0,0 +1,24 @@
+add_input_file_links()
+
+
+add_dumux_test(test_boxrichardsnc test_boxrichardsnc test_boxrichardsnc.cc
+ python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/richardslensbox-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/richardswelltracerbox-00008.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_boxrichardsnc")
+
+add_dumux_test(test_ccrichardsnc test_ccrichardsnc test_ccrichardsnc.cc
+ python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
+ --script fuzzy
+ --files ${CMAKE_SOURCE_DIR}/test/references/richardslenscc-reference.vtu
+ ${CMAKE_CURRENT_BINARY_DIR}/richardswelltracercc-00007.vtu
+ --command "${CMAKE_CURRENT_BINARY_DIR}/test_ccrichardsnc")
+
+#install sources
+install(FILES
+richardswelltracerproblem.hh
+richardswelltracerspatialparams.hh
+test_boxrichardsnc.cc
+test_ccrichardsnc.cc
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/test/porousmediumflow/richardsnc/implicit)