diff --git a/dumux/porousmediumflow/richards/implicit/indices.hh b/dumux/porousmediumflow/richards/implicit/indices.hh
index 32be851c89ffc5fa6d33786d14057ecdb1c42a51..0216062c0c59debbfc10c8bf40b20ae80d7ac116 100644
--- a/dumux/porousmediumflow/richards/implicit/indices.hh
+++ b/dumux/porousmediumflow/richards/implicit/indices.hh
@@ -55,8 +55,8 @@ struct RichardsIndices
     //////////
     // phase indices
     //////////
-    static const int wPhaseIdx = FluidSystem::wPhaseIdx; //!< Index of the wetting phase;
-    static const int nPhaseIdx = FluidSystem::nPhaseIdx; //!< Index of the non-wetting phase;
+    static const int wPhaseIdx = 0; //!< Index of the wetting phase;
+    static const int nPhaseIdx = 1; //!< Index of the non-wetting phase;
 };
 // \}
 
diff --git a/dumux/porousmediumflow/richards/implicit/localresidual.hh b/dumux/porousmediumflow/richards/implicit/localresidual.hh
deleted file mode 100644
index 5659a08f071c772ba1e286759ec405b2d4288d9b..0000000000000000000000000000000000000000
--- a/dumux/porousmediumflow/richards/implicit/localresidual.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 Element-wise calculation of the Jacobian matrix for problems
- *        using the n-phase immiscible fully implicit models.
- */
-#ifndef DUMUX_RICHARDS_LOCAL_RESIDUAL_HH
-#define DUMUX_RICHARDS_LOCAL_RESIDUAL_HH
-
-namespace Dumux
-{
-/*!
- * \ingroup OnePModel
- * \ingroup RichardsLocalResidual
- * \brief Element-wise calculation of the Jacobian matrix for problems
- *        using the n-phase immiscible fully implicit models.
- */
-template<class TypeTag>
-class RichardsLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
-{
-    typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
-
-    using ParentType = typename GET_PROP_TYPE(TypeTag, BaseLocalResidual);
-    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using PrimaryVariables = typename GET_PROP_TYPE(TypeTag, PrimaryVariables);
-    using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
-    using ElementVolumeVariables = typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables);
-    using FluxVariables = typename GET_PROP_TYPE(TypeTag, FluxVariables);
-    using ElementFluxVariablesCache = typename GET_PROP_TYPE(TypeTag, ElementFluxVariablesCache);
-    using SubControlVolume = typename GET_PROP_TYPE(TypeTag, SubControlVolume);
-    using SubControlVolumeFace = typename GET_PROP_TYPE(TypeTag, SubControlVolumeFace);
-    using FVElementGeometry = typename GET_PROP_TYPE(TypeTag, FVElementGeometry);
-    using GridView = typename GET_PROP_TYPE(TypeTag, GridView);
-    using Element = typename GridView::template Codim<0>::Entity;
-    using EnergyLocalResidual = typename GET_PROP_TYPE(TypeTag, EnergyLocalResidual);
-    using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
-    // first index for the mass balance
-    enum { conti0EqIdx = Indices::conti0EqIdx };
-
-    enum {
-        nPhaseIdx = Indices::nPhaseIdx,
-        wPhaseIdx = Indices::wPhaseIdx
-    };
-
-public:
-
-    /*!
-     * \brief Evaluate the rate of change of all conservation
-     *        quantites (e.g. phase mass) within a sub-control
-     *        volume of a finite volume element for the immiscible models.
-     * \param scv The sub control volume
-     * \param volVars The current or previous volVars
-     * \note This function should not include the source and sink terms.
-     * \note The volVars can be different to allow computing
-     *       the implicit euler time derivative here
-     */
-    PrimaryVariables computeStorage(const SubControlVolume& scv,
-                                    const VolumeVariables& volVars) const
-    {
-        // partial time derivative of the phase mass
-        PrimaryVariables storage;
-        storage[conti0EqIdx] = volVars.porosity()
-                               * volVars.density(wPhaseIdx)
-                               * volVars.saturation(wPhaseIdx);
-
-        //! The energy storage in the fluid phase with index phaseIdx
-        EnergyLocalResidual::fluidPhaseStorage(storage, scv, volVars, wPhaseIdx);
-
-        //! The energy storage in the solid matrix
-        EnergyLocalResidual::solidPhaseStorage(storage, scv, volVars);
-
-        return storage;
-    }
-
-
-    /*!
-     * \brief Evaluate the mass flux over a face of a sub control volume
-     * \param scvf The sub control volume face to compute the flux on
-     */
-    PrimaryVariables computeFlux(const Element& element,
-                                 const FVElementGeometry& fvGeometry,
-                                 const ElementVolumeVariables& elemVolVars,
-                                 const SubControlVolumeFace& scvf,
-                                 const ElementFluxVariablesCache& elemFluxVarsCache) const
-    {
-        FluxVariables fluxVars;
-        fluxVars.init(this->problem(), element, fvGeometry, elemVolVars, scvf, elemFluxVarsCache);
-
-        PrimaryVariables flux;
-        // the physical quantities for which we perform upwinding
-        auto upwindTerm = [](const VolumeVariables& volVars)
-                          { return volVars.density(wPhaseIdx)*volVars.mobility(wPhaseIdx); };
-
-        flux[conti0EqIdx] = fluxVars.advectiveFlux(wPhaseIdx, upwindTerm);
-
-        //! Add advective phase energy fluxes. For isothermal model the contribution is zero.
-        EnergyLocalResidual::heatConvectionFlux(flux, fluxVars, wPhaseIdx);
-
-        //! Add diffusive energy fluxes. For isothermal model the contribution is zero.
-        EnergyLocalResidual::heatConductionFlux(flux, fluxVars);
-
-        return flux;
-    }
-
-private:
-    Implementation *asImp_()
-    { return static_cast<Implementation *> (this); }
-
-    const Implementation *asImp_() const
-    { return static_cast<const Implementation *> (this); }
-};
-
-} // end namespace Dumux
-
-#endif
diff --git a/dumux/porousmediumflow/richards/implicit/model.hh b/dumux/porousmediumflow/richards/implicit/model.hh
index b27bc52cfc937f65017737e5eba4ea4768757f56..d71307fc9b43289df4919d2bf9130f07390569bf 100644
--- a/dumux/porousmediumflow/richards/implicit/model.hh
+++ b/dumux/porousmediumflow/richards/implicit/model.hh
@@ -26,7 +26,7 @@
 #define DUMUX_RICHARDS_MODEL_HH
 
 #include <dumux/implicit/model.hh>
-#include <dumux/porousmediumflow/implicit/velocityoutput.hh>
+#include <dumux/porousmediumflow/nonisothermal/implicit/model.hh>
 #include "properties.hh"
 
 namespace Dumux
@@ -108,6 +108,9 @@ class RichardsModel : public GET_PROP_TYPE(TypeTag, BaseModel)
     using Problem = typename GET_PROP_TYPE(TypeTag, Problem);
     using VolumeVariables = typename GET_PROP_TYPE(TypeTag, VolumeVariables);
     using Indices = typename GET_PROP_TYPE(TypeTag, Indices);
+
+    using NonIsothermalModel = Dumux::NonIsothermalModel<TypeTag>;
+
     enum {
         nPhaseIdx = Indices::nPhaseIdx,
         wPhaseIdx = Indices::wPhaseIdx
@@ -132,15 +135,15 @@ public:
         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("rhoW", [](const VolumeVariables& v){ return v.density(wPhaseIdx); });
-        vtkOutputModule.addSecondaryVariable("rhoN", [](const VolumeVariables& v){ return v.density(nPhaseIdx); });
-        vtkOutputModule.addSecondaryVariable("mobW", [](const VolumeVariables& v){ return v.mobility(wPhaseIdx); });
-        vtkOutputModule.addSecondaryVariable("mobN", [](const VolumeVariables& v){ return v.mobility(nPhaseIdx); });
-        vtkOutputModule.addSecondaryVariable("temperature", [](const VolumeVariables& v){ return v.temperature(); });
+        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("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); });
+
+        NonIsothermalModel::maybeAddTemperature(vtkOutputModule);
     }
 };
 
diff --git a/dumux/porousmediumflow/richards/implicit/propertydefaults.hh b/dumux/porousmediumflow/richards/implicit/propertydefaults.hh
index c1948c41c26381dff16c69eb2177f54cf9fc8b11..859ca527893b3e073554ef2813f1270c3ed22d9e 100644
--- a/dumux/porousmediumflow/richards/implicit/propertydefaults.hh
+++ b/dumux/porousmediumflow/richards/implicit/propertydefaults.hh
@@ -32,14 +32,15 @@
 #include "indices.hh"
 #include "volumevariables.hh"
 #include "properties.hh"
-#include "localresidual.hh"
 #include "newtoncontroller.hh"
 
+#include <dumux/porousmediumflow/immiscible/localresidual.hh>
 #include <dumux/porousmediumflow/nonisothermal/implicit/propertydefaults.hh>
-#include <dumux/material/components/nullcomponent.hh>
 #include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
-#include <dumux/material/fluidsystems/2pimmiscible.hh>
 #include <dumux/material/spatialparams/implicit.hh>
+#include <dumux/material/components/simpleh2o.hh>
+#include <dumux/material/fluidsystems/liquidphase.hh>
+#include <dumux/material/fluidstates/immiscible.hh>
 
 namespace Dumux
 {
@@ -51,11 +52,15 @@ namespace Properties {
 //////////////////////////////////////////////////////////////////
 //! Number of equations required by the model
 SET_INT_PROP(Richards, NumEq, 1);
+
 //! Number of fluid phases considered
-SET_INT_PROP(Richards, NumPhases, 2);
+//! Although the number of phases is two for the Richards model the
+//! non-wetting phase is not balanced and thus we only have one
+//! phase with a balance equation
+SET_INT_PROP(Richards, NumPhases, 1);
 
 //! The local residual operator
-SET_TYPE_PROP(Richards, LocalResidual, RichardsLocalResidual<TypeTag>);
+SET_TYPE_PROP(Richards, LocalResidual, ImmiscibleLocalResidual<TypeTag>);
 
 //! The global model used
 SET_TYPE_PROP(Richards, Model, RichardsModel<TypeTag>);
@@ -89,59 +94,15 @@ SET_TYPE_PROP(Richards, SpatialParams, ImplicitSpatialParams<TypeTag>);
  */
 SET_TYPE_PROP(Richards, MaterialLawParams, typename GET_PROP_TYPE(TypeTag, MaterialLaw)::Params);
 
-/*!
- * \brief The wetting phase used.
- *
- * By default we use the null-phase, i.e. this has to be defined by
- * the problem for the program to work. Please be aware that you
- * should be careful to use the Richards model in conjunction with
- * liquid non-wetting phases. This is only meaningful if the viscosity
- * of the liquid phase is _much_ lower than the viscosity of the
- * wetting phase.
- */
-
-SET_PROP(Richards, WettingPhase)
-{
-private:
-    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-public:
-    using type = FluidSystems::LiquidPhase<Scalar, NullComponent<Scalar>>;
-};
-
-/*!
- * \brief The wetting phase used.
- *
- * By default we use the null-phase, i.e. this has to be defined by
- * the problem for the program to work. This doed not need to be
- * specified by the problem for the Richards model to work because the
- * Richards model does not conserve the non-wetting phase.
- */
-SET_PROP(Richards, NonwettingPhase)
-{
-private:
-    using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-public:
-    using type = FluidSystems::GasPhase<Scalar, NullComponent<Scalar>>;
-};
-
 /*!
  *\brief The fluid system used by the model.
  *
- * By default this uses the immiscible twophase fluid system. The
- * actual fluids used are specified using in the problem definition by
- * the WettingPhase and NonwettingPhase properties. Be aware that
- * using different fluid systems in conjunction with the Richards
- * model only makes very limited sense.
+ * By default this uses the liquid phase fluid system with simple H2O.
  */
 SET_PROP(Richards, FluidSystem)
 {
-private:
     using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar);
-    using WettingPhase = typename GET_PROP_TYPE(TypeTag, WettingPhase);
-    using NonwettingPhase = typename GET_PROP_TYPE(TypeTag, NonwettingPhase);
-
-public:
-    using type = FluidSystems::TwoPImmiscible<Scalar, WettingPhase, NonwettingPhase>;
+    using type = FluidSystems::LiquidPhase<Scalar, SimpleH2O<Scalar>>;
 };
 
 /*!
@@ -190,7 +151,7 @@ SET_TYPE_PROP(RichardsNI, IsothermalModel, RichardsModel<TypeTag>);
 SET_TYPE_PROP(RichardsNI, IsothermalVolumeVariables, RichardsVolumeVariables<TypeTag>);
 
 //set isothermal LocalResidual
-SET_TYPE_PROP(RichardsNI, IsothermalLocalResidual, RichardsLocalResidual<TypeTag>);
+SET_TYPE_PROP(RichardsNI, IsothermalLocalResidual, ImmiscibleLocalResidual<TypeTag>);
 
 //set isothermal Indices
 SET_TYPE_PROP(RichardsNI, IsothermalIndices, RichardsIndices<TypeTag>);
diff --git a/dumux/porousmediumflow/richards/implicit/volumevariables.hh b/dumux/porousmediumflow/richards/implicit/volumevariables.hh
index 2c5eca5feabc627b1c1384db869f67846a515bbe..ddf3d0372e880fd15fdf923c4c749c3c98ec5d3a 100644
--- a/dumux/porousmediumflow/richards/implicit/volumevariables.hh
+++ b/dumux/porousmediumflow/richards/implicit/volumevariables.hh
@@ -26,7 +26,6 @@
 
 #include "properties.hh"
 
-#include <dumux/material/fluidstates/immiscible.hh>
 #include <dumux/discretization/volumevariables.hh>
 
 namespace Dumux
@@ -79,8 +78,6 @@ public:
                 const Element &element,
                 const SubControlVolume& scv)
     {
-        assert(!FluidSystem::isLiquid(nPhaseIdx));
-
         ParentType::update(elemSol, problem, element, scv);
 
         completeFluidState(elemSol, problem, element, scv, fluidState_);
@@ -93,9 +90,6 @@ public:
         pnRef_ = problem.nonWettingReferencePressure();
         porosity_ = problem.spatialParams().porosity(element, scv, elemSol);
         permeability_ = problem.spatialParams().permeability(element, scv, elemSol);
-
-        // energy related quantities not belonging to the fluid state
-        asImp_().updateEnergy_(elemSol, problem, element, scv);
     }
 
     /*!
@@ -107,7 +101,7 @@ public:
                                    const SubControlVolume& scv,
                                    FluidState& fluidState)
     {
-        Scalar t = Implementation::temperature_(elemSol, problem, element, scv);
+        Scalar t = ParentType::temperature(elemSol, problem, element, scv);
         fluidState.setTemperature(t);
 
         const auto& materialParams = problem.spatialParams().materialLawParams(element, scv, elemSol);
@@ -115,26 +109,23 @@ public:
 
         const Scalar minPc = MaterialLaw::pc(materialParams, 1.0);
         fluidState.setPressure(wPhaseIdx, priVars[pressureIdx]);
-        fluidState.setPressure(nPhaseIdx, std::max(problem.nonWettingReferencePressure(), priVars[pressureIdx] + minPc));
+
+        using std::max;
+        pn_ = max(problem.nonWettingReferencePressure(), priVars[pressureIdx] + minPc);
 
         // saturations
-        const Scalar sw = MaterialLaw::sw(materialParams, fluidState.pressure(nPhaseIdx) - fluidState.pressure(wPhaseIdx));
+        const Scalar sw = MaterialLaw::sw(materialParams, pn_ - fluidState.pressure(wPhaseIdx));
         fluidState.setSaturation(wPhaseIdx, sw);
-        fluidState.setSaturation(nPhaseIdx, 1 - sw);
 
         // density and viscosity
         typename FluidSystem::ParameterCache paramCache;
         paramCache.updateAll(fluidState);
         fluidState.setDensity(wPhaseIdx, FluidSystem::density(fluidState, paramCache, wPhaseIdx));
-        fluidState.setDensity(nPhaseIdx, 1e-10);
 
         fluidState.setViscosity(wPhaseIdx, FluidSystem::viscosity(fluidState, paramCache, wPhaseIdx));
-        fluidState.setViscosity(nPhaseIdx, 1e-10);
 
         // compute and set the enthalpy
-        fluidState.setEnthalpy(wPhaseIdx, Implementation::enthalpy_(fluidState, paramCache, wPhaseIdx));
-        fluidState.setEnthalpy(nPhaseIdx, Implementation::enthalpy_(fluidState, paramCache, nPhaseIdx));
-
+        fluidState.setEnthalpy(wPhaseIdx, Implementation::enthalpy(fluidState, paramCache, wPhaseIdx));
     }
 
     /*!
@@ -144,6 +135,12 @@ public:
     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.
      *
@@ -170,7 +167,7 @@ public:
      * \param phaseIdx The index of the fluid phase
      */
     Scalar saturation(const int phaseIdx) const
-    { return fluidState_.saturation(phaseIdx); }
+    { 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
@@ -179,7 +176,7 @@ public:
      * \param phaseIdx The index of the fluid phase
      */
     Scalar density(const int phaseIdx) const
-    { return fluidState_.density(phaseIdx); }
+    { return phaseIdx == wPhaseIdx ? fluidState_.density(phaseIdx) : 0.0; }
 
     /*!
      * \brief Returns the effective pressure \f$\mathrm{[Pa]}\f$ of a given phase within
@@ -193,17 +190,7 @@ public:
      * \param phaseIdx The index of the fluid phase
      */
     Scalar pressure(const int phaseIdx) const
-    { return fluidState_.pressure(phaseIdx); }
-
-    /*!
-     * \brief Returns average temperature \f$\mathrm{[K]}\f$ inside the control volume.
-     *
-     * Note that we assume thermodynamic equilibrium, i.e. the
-     * temperature of the rock matrix and of all fluid phases are
-     * identical.
-     */
-    Scalar temperature() const
-    { return fluidState_.temperature(); }
+    { return phaseIdx == wPhaseIdx ? fluidState_.pressure(phaseIdx) : pn_; }
 
     /*!
      * \brief Returns the effective mobility \f$\mathrm{[1/(Pa*s)]}\f$ of a given phase within
@@ -219,6 +206,16 @@ public:
     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.
@@ -226,11 +223,7 @@ public:
      * \param phaseIdx The index of the fluid phase
      */
     Scalar relativePermeability(const int phaseIdx) const
-    {
-        if (phaseIdx == wPhaseIdx)
-            return relativePermeabilityWetting_;
-        return 1;
-    }
+    { return phaseIdx == wPhaseIdx ? relativePermeabilityWetting_ : 1.0; }
 
     /*!
      * \brief Returns the effective capillary pressure \f$\mathrm{[Pa]}\f$ within the
@@ -241,9 +234,7 @@ public:
      * \f[ p_c = p_n - p_w \f]
      */
     Scalar capillaryPressure() const
-    {
-        return fluidState_.pressure(nPhaseIdx) - fluidState_.pressure(wPhaseIdx);
-    }
+    { return pn_ - fluidState_.pressure(wPhaseIdx); }
 
     /*!
      * \brief Returns the pressureHead \f$\mathrm{[cm]}\f$ of a given phase within
@@ -260,9 +251,7 @@ public:
      *       or the gravity different
      */
     Scalar pressureHead(const int phaseIdx) const
-    {
-        return 100.0 *(fluidState_.pressure(phaseIdx) - pnRef_)/fluidState_.density(phaseIdx)/9.81;
-    }
+    { return 100.0 *(pressure(phaseIdx) - pnRef_)/density(phaseIdx)/9.81; }
 
     /*!
      * \brief Returns the water content
@@ -276,41 +265,15 @@ public:
      *       manually do a conversion.
      */
     Scalar waterContent(const int phaseIdx) const
-    {
-        return fluidState_.saturation(phaseIdx)* porosity_;
-    }
+    { return saturation(phaseIdx) * porosity_; }
 
 protected:
-    static Scalar temperature_(const ElementSolutionVector &elemSol,
-                               const Problem& problem,
-                               const Element &element,
-                               const SubControlVolume &scv)
-    {
-        return problem.temperatureAtPos(scv.dofPosition());
-    }
-
-    template<class ParameterCache>
-    static Scalar enthalpy_(const FluidState& fluidState,
-                            const ParameterCache& paramCache,
-                            const int phaseIdx)
-    {
-        return 0;
-    }
-
-    /*!
-     * \brief Called by update() to compute the energy related quantities.
-     */
-    void updateEnergy_(const ElementSolutionVector &elemSol,
-                       const Problem &problem,
-                       const Element &element,
-                       const SubControlVolume& scv)
-    {}
-
     FluidState fluidState_;
     Scalar relativePermeabilityWetting_;
     Scalar porosity_;
     PermeabilityType permeability_;
     Scalar pnRef_;
+    static Scalar pn_;
 
 private:
     Implementation &asImp_()
@@ -320,6 +283,9 @@ private:
     { return *static_cast<const Implementation*>(this); }
 };
 
+template <class TypeTag>
+typename GET_PROP_TYPE(TypeTag, Scalar) RichardsVolumeVariables<TypeTag>::pn_;
+
 } // end namespace Dumux
 
 #endif
diff --git a/test/porousmediumflow/richards/implicit/CMakeLists.txt b/test/porousmediumflow/richards/implicit/CMakeLists.txt
index 96adf2a40495c9f373d9a62c8769d8cca6c36d7b..0726f6c2a829274589ac21a57f65a73bb55abdfd 100644
--- a/test/porousmediumflow/richards/implicit/CMakeLists.txt
+++ b/test/porousmediumflow/richards/implicit/CMakeLists.txt
@@ -5,7 +5,7 @@ if(MPI_FOUND)
                  python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
                    --script fuzzy
                    --files ${CMAKE_SOURCE_DIR}/test/references/richardslensbox-reference-parallel.vtu
-                           ${CMAKE_CURRENT_BINARY_DIR}/s0002-p0000-richardslensbox-00007.vtu
+                           ${CMAKE_CURRENT_BINARY_DIR}/s0002-p0000-richardslensbox-00008.vtu
                    --command "${MPIEXEC} -np 2 ${CMAKE_CURRENT_BINARY_DIR}/test_boxrichards")
 else()
 # isothermal tests
@@ -13,7 +13,7 @@ else()
                  python ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py
                    --script fuzzy
                    --files ${CMAKE_SOURCE_DIR}/test/references/richardslensbox-reference.vtu
-                           ${CMAKE_CURRENT_BINARY_DIR}/richardslensbox-00007.vtu
+                           ${CMAKE_CURRENT_BINARY_DIR}/richardslensbox-00008.vtu
                    --command "${CMAKE_CURRENT_BINARY_DIR}/test_boxrichards")
 endif()
 
diff --git a/test/porousmediumflow/richards/implicit/richardslensspatialparams.hh b/test/porousmediumflow/richards/implicit/richardslensspatialparams.hh
index 6c589ceac44bdae7917d20d8087a914163cc311b..19a6914b5a36b87a117a0dc3b374bc674f1baeee 100644
--- a/test/porousmediumflow/richards/implicit/richardslensspatialparams.hh
+++ b/test/porousmediumflow/richards/implicit/richardslensspatialparams.hh
@@ -95,7 +95,7 @@ public:
     {
 
         lensLowerLeft_ = {1.0, 2.0};
-        lensUpperRight_ = {4.0, 3.0};
+        lensUpperRight_ = {4.1, 3.1};
 
         // residual saturations
         lensMaterialParams_.setSwr(0.18);
diff --git a/test/references/richardslensbox-reference-parallel.vtu b/test/references/richardslensbox-reference-parallel.vtu
index acca99ab726d1138e223c8b4e9deff0097b47aa7..d64d207ecaab3bec5a74bec05290d0630128f372 100644
--- a/test/references/richardslensbox-reference-parallel.vtu
+++ b/test/references/richardslensbox-reference-parallel.vtu
@@ -2,28 +2,7 @@
 <VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
   <UnstructuredGrid>
     <Piece NumberOfCells="192" NumberOfPoints="221">
-      <PointData Scalars="Sn">
-        <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii">
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 0.999972 0.999959
-          0.999972 1 1 1 1 1 1 1 1 1 0.997392 0.909005
-          0.904424 0.909005 0.997392 1 1 1 1 1 1 1 0.998941 0.771659
-          0.385726 0.365066 0.385726 0.771659 1 1 1 1 1 1 0.999996 0.918448
-          0.287481 0.0673969 0.0636792 0.0673969 0.287481
-        </DataArray>
+      <PointData Scalars="Sw">
         <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii">
           9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
           9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
@@ -45,26 +24,26 @@
           0.614274 0.634934 0.614274 0.228341 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 4.31542e-06 0.0815522
           0.712519 0.932603 0.936321 0.932603 0.712519
         </DataArray>
-        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000
+        <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii">
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 0.999972 0.999959
+          0.999972 1 1 1 1 1 1 1 1 1 0.997392 0.909005
+          0.904424 0.909005 0.997392 1 1 1 1 1 1 1 0.998941 0.771659
+          0.385726 0.365066 0.385726 0.771659 1 1 1 1 1 1 0.999996 0.918448
+          0.287481 0.0673969 0.0636792 0.0673969 0.287481
         </DataArray>
         <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
           97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4
@@ -87,6 +66,27 @@
           99733.4 99738.8 99733.4 99586.6 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.5 99323.6
           99759.2 99835.1 99837.3 99835.1 99759.2
         </DataArray>
+        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000
+        </DataArray>
         <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
           2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62
           2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62
@@ -108,7 +108,7 @@
           266.621 261.187 266.621 413.416 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.51 676.417
           240.748 164.873 162.729 164.873 240.748
         </DataArray>
-        <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="density" NumberOfComponents="1" format="ascii">
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
@@ -129,28 +129,7 @@
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000
         </DataArray>
-        <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii">
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10
-        </DataArray>
-        <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="mobility" NumberOfComponents="1" format="ascii">
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
@@ -171,26 +150,26 @@
           145.944 164.193 145.944 3.93436 0 0 0 0 0 0 0 0.0198338
           248.432 732.351 752.186 732.351 248.432
         </DataArray>
-        <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii">
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10
+        <DataArray type="Float32" Name="kr" NumberOfComponents="1" format="ascii">
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 1.57172e-05
+          2.3053e-05 1.57172e-05 0 0 0 0 0 0 0 0 0 0.00287218
+          0.120485 0.135536 0.120485 0.00287218 0 0 0 0 0 0 0 1.05648e-05
+          0.237214 0.71697 0.737982 0.71697 0.237214
         </DataArray>
         <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
           0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
@@ -213,27 +192,6 @@
           0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
           0.4 0.4 0.4 0.4 0.4
         </DataArray>
-        <DataArray type="Float32" Name="temperature" NumberOfComponents="1" format="ascii">
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15
-        </DataArray>
         <DataArray type="Float32" Name="pressure head" NumberOfComponents="1" format="ascii">
           -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861
           -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861
diff --git a/test/references/richardslensbox-reference.vtu b/test/references/richardslensbox-reference.vtu
index 9dff729c81486ee809123f190e76fe018b50afbc..2d03ab8ffe8f5c0f22694fbcdb2f62a3d1df9a56 100644
--- a/test/references/richardslensbox-reference.vtu
+++ b/test/references/richardslensbox-reference.vtu
@@ -2,45 +2,7 @@
 <VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
   <UnstructuredGrid>
     <Piece NumberOfCells="384" NumberOfPoints="425">
-      <PointData Scalars="Sn">
-        <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii">
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 0.999976 0.999964
-          0.999976 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 0.997534 0.909917
-          0.905361 0.909917 0.997534 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 0.999 0.771847
-          0.384776 0.36421 0.384776 0.771847 0.999 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 0.999996 0.918527
-          0.286966 0.0671314 0.0634913 0.0671314 0.286966 0.918527 0.999996 1 1 1 1 1
-          1 1 1 1 1
-        </DataArray>
+      <PointData Scalars="Sw">
         <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii">
           9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
           9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
@@ -79,43 +41,43 @@
           0.713034 0.932869 0.936509 0.932869 0.713034 0.0814728 3.93768e-06 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
           9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
         </DataArray>
-        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000
+        <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii">
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 0.999976 0.999964
+          0.999976 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 0.997534 0.909917
+          0.905361 0.909917 0.997534 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 0.999 0.771847
+          0.384776 0.36421 0.384776 0.771847 0.999 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 0.999996 0.918527
+          0.286966 0.0671314 0.0634913 0.0671314 0.286966 0.918527 0.999996 1 1 1 1 1
+          1 1 1 1 1
         </DataArray>
         <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
           97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4
@@ -155,6 +117,44 @@
           99759.4 99835.3 99837.4 99835.3 99759.4 99323.1 97470.5 97470.4 97470.4 97470.4 97470.4 97470.4
           97470.4 97470.4 97470.4 97470.4 97470.4
         </DataArray>
+        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000
+        </DataArray>
         <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
           2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62
           2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62
@@ -193,7 +193,7 @@
           240.611 164.722 162.618 164.722 240.611 676.884 2529.52 2529.62 2529.62 2529.62 2529.62 2529.62
           2529.62 2529.62 2529.62 2529.62 2529.62
         </DataArray>
-        <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="density" NumberOfComponents="1" format="ascii">
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
@@ -231,45 +231,7 @@
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000
         </DataArray>
-        <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii">
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10
-        </DataArray>
-        <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="mobility" NumberOfComponents="1" format="ascii">
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
@@ -307,44 +269,6 @@
           249.082 733.751 753.199 733.751 249.082 0.0196823 0 0 0 0 0 0
           0 0 0 0 0
         </DataArray>
-        <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii">
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10
-        </DataArray>
         <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
           0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
           0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
@@ -383,44 +307,6 @@
           0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
           0.4 0.4 0.4 0.4 0.4
         </DataArray>
-        <DataArray type="Float32" Name="temperature" NumberOfComponents="1" format="ascii">
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15
-        </DataArray>
         <DataArray type="Float32" Name="pressure head" NumberOfComponents="1" format="ascii">
           -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861
           -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861
diff --git a/test/references/richardslenscc-reference.vtu b/test/references/richardslenscc-reference.vtu
index cca30ff2454307aafdd825ff8fe361a9f4521d58..fde4bda5f6548856002d86db1cbc1b9fe5812bd6 100644
--- a/test/references/richardslenscc-reference.vtu
+++ b/test/references/richardslenscc-reference.vtu
@@ -2,41 +2,7 @@
 <VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
   <UnstructuredGrid>
     <Piece NumberOfCells="384" NumberOfPoints="425">
-      <CellData Scalars="Sn">
-        <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii">
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 1 0.992854 0.988373 0.988373 0.992854
-          1 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 1 0.996667 0.703082 0.660102 0.660102 0.703082
-          0.996667 1 1 1 1 1 1 1 1 1 1 1
-          1 1 1 1 1 1 0.999997 0.920927 0.185263 0.152734 0.152734 0.185263
-          0.920927 0.999997 1 1 1 1 1 1 1 1 1 1
-        </DataArray>
+      <CellData Scalars="Sw">
         <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii">
           9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
           9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
@@ -71,39 +37,39 @@
           9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 2.85357e-06 0.0790732 0.814737 0.847266 0.847266 0.814737
           0.0790732 2.85357e-06 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17 9.02056e-17
         </DataArray>
-        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
-          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+        <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii">
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 0.992854 0.988373 0.988373 0.992854
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 0.996667 0.703082 0.660102 0.660102 0.703082
+          0.996667 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 0.999997 0.920927 0.185263 0.152734 0.152734 0.185263
+          0.920927 0.999997 1 1 1 1 1 1 1 1 1 1
         </DataArray>
         <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
           97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4
@@ -139,6 +105,40 @@
           97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.5 99308.3 99788.3 99798.9 99798.9 99788.3
           99308.3 97470.5 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4 97470.4
         </DataArray>
+        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+          100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
+        </DataArray>
         <DataArray type="Float32" Name="pc" NumberOfComponents="1" format="ascii">
           2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62
           2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62
@@ -173,7 +173,7 @@
           2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.54 691.737 211.699 201.127 201.127 211.699
           691.737 2529.54 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62 2529.62
         </DataArray>
-        <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="density" NumberOfComponents="1" format="ascii">
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
@@ -207,41 +207,7 @@
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
           1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000
         </DataArray>
-        <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii">
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-          1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10 1e-10
-        </DataArray>
-        <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="mobility" NumberOfComponents="1" format="ascii">
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
@@ -275,39 +241,39 @@
           0 0 0 0 0 0 0 0.0154609 407.511 473.279 473.279 407.511
           0.0154609 0 0 0 0 0 0 0 0 0 0 0
         </DataArray>
-        <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii">
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
-          1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10 1e+10
+        <DataArray type="Float32" Name="kr" NumberOfComponents="1" format="ascii">
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0.0107357 0.0176661 0.0176661 0.0107357
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 1.54609e-05 0.407511 0.473279 0.473279 0.407511
+          1.54609e-05 0 0 0 0 0 0 0 0 0 0 0
         </DataArray>
         <DataArray type="Float32" Name="porosity" NumberOfComponents="1" format="ascii">
           0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
@@ -343,40 +309,6 @@
           0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
           0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4
         </DataArray>
-        <DataArray type="Float32" Name="temperature" NumberOfComponents="1" format="ascii">
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-          283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15 283.15
-        </DataArray>
         <DataArray type="Float32" Name="pressure head" NumberOfComponents="1" format="ascii">
           -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861
           -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861 -25.7861