diff --git a/dumux/implicit/richards/richardsfluxvariables.hh b/dumux/implicit/richards/richardsfluxvariables.hh
new file mode 100644
index 0000000000000000000000000000000000000000..131607302c396f790ff86989278dc1d2f344391e
--- /dev/null
+++ b/dumux/implicit/richards/richardsfluxvariables.hh
@@ -0,0 +1,108 @@
+// -*- 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 Overwrites the volumeFlux() function from the ImplicitDarcyFluxVariables.
+ *
+ */
+#ifndef DUMUX_RICHARDS_FLUX_VARIABLES_HH
+#define DUMUX_RICHARDS_FLUX_VARIABLES_HH
+
+#include <dumux/common/math.hh>
+#include <dumux/common/parameters.hh>
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
+#include "richardsproperties.hh"
+
+namespace Dumux
+{
+    
+
+/*!
+ * \ingroup RichardsModel
+ * \ingroup ImplicitFluxVariables
+ * \brief Overwrites the volumeFlux() function from the ImplicitDarcyFluxVariables.
+ */
+template <class TypeTag>
+class RichardsFluxVariables : public ImplicitDarcyFluxVariables<TypeTag>
+{
+  
+    typedef ImplicitDarcyFluxVariables<TypeTag> ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;   
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Problem) Problem;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+    typedef typename GridView::template Codim<0>::Entity Element;
+
+    enum { dim = GridView::dimension} ;
+    enum { dimWorld = GridView::dimensionworld} ;
+    enum { numPhases = GET_PROP_VALUE(TypeTag, NumPhases)} ;
+    enum { nPhaseIdx = Indices::nPhaseIdx} ;
+
+
+
+public:
+    /*
+     * \brief The constructor
+     *
+     * \param problem The problem
+     * \param element The finite element
+     * \param fvGeometry The finite-volume geometry
+     * \param faceIdx The local index of the SCV (sub-control-volume) face
+     * \param elemVolVars The volume variables of the current element
+     * \param onBoundary A boolean variable to specify whether the flux variables
+     * are calculated for interior SCV faces or boundary faces, default=false
+     */
+    RichardsFluxVariables(const Problem &problem,
+                 const Element &element,
+                 const FVElementGeometry &fvGeometry,
+                 const int fIdx,
+                 const ElementVolumeVariables &elemVolVars,
+                 const bool onBoundary = false)
+    : ParentType(problem, element, fvGeometry, fIdx, elemVolVars, onBoundary)
+    {    }
+
+public:
+    /*!
+     * \brief Return the volumetric flux over a face of a given phase.
+     *
+     *        This is the calculated velocity multiplied by the unit normal
+     *        and the area of the face.
+     *        face().normal
+     *        has already the magnitude of the area. 
+     *        For the Richards model the velocity of the non-wetting phase
+     *        is set to zero.
+     *
+     * \param phaseIdx index of the phase
+     */
+    Scalar volumeFlux(const unsigned int phaseIdx) const
+    { 
+      if(phaseIdx == nPhaseIdx)
+          return 0.;
+      else
+          return ParentType::volumeFlux(phaseIdx);
+    }
+    
+ 
+};
+
+} // end namespace
+
+#endif
diff --git a/dumux/implicit/richards/richardslocalresidual.hh b/dumux/implicit/richards/richardslocalresidual.hh
index c57a514cb700c3716db3c1b94a7aa7ef4e143014..14918e5859ad4e22d2adedb61889f01032b0ad7e 100644
--- a/dumux/implicit/richards/richardslocalresidual.hh
+++ b/dumux/implicit/richards/richardslocalresidual.hh
@@ -36,6 +36,7 @@ namespace Dumux
 template<class TypeTag>
 class RichardsLocalResidual : public GET_PROP_TYPE(TypeTag, BaseLocalResidual)
 {
+    typedef typename GET_PROP_TYPE(TypeTag, LocalResidual) Implementation;
     typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
     typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
     typedef typename GET_PROP_TYPE(TypeTag, FluxVariables) FluxVariables;
@@ -113,12 +114,27 @@ public:
                                this->curVolVars_(),
                                onBoundary);
 
+        flux = 0;
+        asImp_()->computeAdvectiveFlux(flux, fluxVars);
+        asImp_()->computeDiffusiveFlux(flux, fluxVars);
+    }
+
+    /*!
+     * \brief Evaluates the advective mass flux of all components over
+     *        a face of a sub-control volume.
+     *
+     * \param flux The advective flux over the sub-control-volume face for each component
+     * \param fluxVars The flux variables at the current SCV
+     */
+
+    void computeAdvectiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
+    {
         // data attached to upstream and the downstream vertices
         // of the current phase
         const VolumeVariables &up = this->curVolVars_(fluxVars.upstreamIdx(wPhaseIdx));
         const VolumeVariables &dn = this->curVolVars_(fluxVars.downstreamIdx(wPhaseIdx));
 
-        flux[contiEqIdx] =
+        flux[contiEqIdx] +=
             fluxVars.volumeFlux(wPhaseIdx)
             *
             ((    massUpwindWeight_)*up.density(wPhaseIdx)
@@ -126,6 +142,24 @@ public:
              (1 - massUpwindWeight_)*dn.density(wPhaseIdx));
     }
 
+    /*!
+     * \brief Adds the diffusive flux to the flux vector over
+     *        the face of a sub-control volume.
+     *
+     * \param flux The diffusive flux over the sub-control-volume face for each phase
+     * \param fluxVars The flux variables at the current SCV
+     *
+     * This function doesn't do anything but may be used by the
+     * non-isothermal three-phase models to calculate diffusive heat
+     * fluxes
+     */
+    void computeDiffusiveFlux(PrimaryVariables &flux, const FluxVariables &fluxVars) const
+    {
+        // diffusive fluxes
+        flux += 0.0;
+    }
+
+
     /*!
      * \brief Calculate the source term of the equation
      *
@@ -143,8 +177,21 @@ public:
                                      this->curVolVars_());
     }
 
+protected:
+    Implementation *asImp_()
+    {
+        return static_cast<Implementation *> (this);
+    }
+
+    const Implementation *asImp_() const
+    {
+        return static_cast<const Implementation *> (this);
+    }
+
 private:
     Scalar massUpwindWeight_;
+
+
 };
 
 }
diff --git a/dumux/implicit/richards/richardsproperties.hh b/dumux/implicit/richards/richardsproperties.hh
index 3fb082ebc6814aaf3e9c0cdc27e34999e9a3a6c2..c390b7af808997693c1a253dabe067ac43e9a580 100644
--- a/dumux/implicit/richards/richardsproperties.hh
+++ b/dumux/implicit/richards/richardsproperties.hh
@@ -30,6 +30,7 @@
 
 #include <dumux/implicit/box/boxproperties.hh>
 #include <dumux/implicit/cellcentered/ccproperties.hh>
+#include <dumux/implicit/nonisothermal/niproperties.hh>
 
 namespace Dumux
 {
@@ -43,11 +44,16 @@ namespace Properties {
 // Type tags
 //////////////////////////////////////////////////////////////////
 
-//! The type tags for the implicit Richards problems
+//! The type tags for the implicit isothermal one-phase two-component problems
 NEW_TYPE_TAG(Richards);
 NEW_TYPE_TAG(BoxRichards, INHERITS_FROM(BoxModel, Richards));
 NEW_TYPE_TAG(CCRichards, INHERITS_FROM(CCModel, Richards));
 
+//! The type tags for the corresponding non-isothermal problems
+NEW_TYPE_TAG(RichardsNI, INHERITS_FROM(Richards, NonIsothermal));
+NEW_TYPE_TAG(BoxRichardsNI, INHERITS_FROM(BoxModel, RichardsNI));
+NEW_TYPE_TAG(CCRichardsNI, INHERITS_FROM(CCModel, RichardsNI));
+
 //////////////////////////////////////////////////////////////////
 // Property tags
 //////////////////////////////////////////////////////////////////
diff --git a/dumux/implicit/richards/richardspropertydefaults.hh b/dumux/implicit/richards/richardspropertydefaults.hh
index 941238a19ab43f481be6909bb4143f47dff65206..aff4cd2956a220ae5d77ec2665f0fb03a2585d77 100644
--- a/dumux/implicit/richards/richardspropertydefaults.hh
+++ b/dumux/implicit/richards/richardspropertydefaults.hh
@@ -28,6 +28,7 @@
 #ifndef DUMUX_RICHARDS_PROPERTY_DEFAULTS_HH
 #define DUMUX_RICHARDS_PROPERTY_DEFAULTS_HH
 
+#include "richardsfluxvariables.hh"
 #include "richardsmodel.hh"
 #include "richardsproblem.hh"
 #include "richardsindices.hh"
@@ -36,10 +37,12 @@
 #include "richardsnewtoncontroller.hh"
 #include "richardslocalresidual.hh"
 
+#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
+#include <dumux/implicit/nonisothermal/nipropertydefaults.hh>
 #include <dumux/material/components/nullcomponent.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
 #include <dumux/material/fluidsystems/2pimmisciblefluidsystem.hh>
 #include <dumux/material/spatialparams/implicitspatialparams.hh>
-#include <dumux/implicit/common/implicitdarcyfluxvariables.hh>
 
 namespace Dumux
 {
@@ -50,23 +53,23 @@ namespace Properties {
 // Properties values
 //////////////////////////////////////////////////////////////////
 //! Number of equations required by the model
-SET_INT_PROP(Richards, NumEq, 1);
+SET_INT_PROP(Richards, NumEq, GET_PROP_VALUE(TypeTag, IsothermalNumEq));
 //! Number of fluid phases considered
 SET_INT_PROP(Richards, NumPhases, 2);
 
 //! The local residual operator
 SET_TYPE_PROP(Richards,
               LocalResidual,
-              RichardsLocalResidual<TypeTag>);
+              typename GET_PROP_TYPE(TypeTag, IsothermalLocalResidual));
 
 //! The global model used
-SET_TYPE_PROP(Richards, Model, RichardsModel<TypeTag>);
+SET_TYPE_PROP(Richards, Model, typename GET_PROP_TYPE(TypeTag, IsothermalModel));
 
 //! The class for the volume averaged quantities
-SET_TYPE_PROP(Richards, VolumeVariables, RichardsVolumeVariables<TypeTag>);
+SET_TYPE_PROP(Richards, VolumeVariables, typename GET_PROP_TYPE(TypeTag, IsothermalVolumeVariables));
 
 //! The class for the quantities required for the flux calculation
-SET_TYPE_PROP(Richards, FluxVariables, ImplicitDarcyFluxVariables<TypeTag>);
+SET_TYPE_PROP(Richards, FluxVariables, typename GET_PROP_TYPE(TypeTag, IsothermalFluxVariables));
 
 //! The class of the newton controller
 SET_TYPE_PROP(Richards, NewtonController, RichardsNewtonController<TypeTag>);
@@ -78,7 +81,7 @@ SET_SCALAR_PROP(Richards, ImplicitMassUpwindWeight, 1.0);
 SET_SCALAR_PROP(Richards, ImplicitMobilityUpwindWeight, 1.0);
 
 //! The class with all index definitions for the model
-SET_TYPE_PROP(Richards, Indices, RichardsIndices<TypeTag>);
+SET_TYPE_PROP(Richards, Indices, typename GET_PROP_TYPE(TypeTag, IsothermalIndices));
 
 //! The spatial parameters to be employed. 
 //! Use ImplicitSpatialParams by default.
@@ -164,6 +167,41 @@ SET_BOOL_PROP(Richards, ProblemEnableGravity, true);
 //        (Nield, Bejan, Convection in porous media, 2006, p. 10)
 SET_SCALAR_PROP(BoxModel, SpatialParamsForchCoeff, 0.55);
 
+//! Somerton is used as default model to compute the effective thermal heat conductivity
+SET_PROP(NonIsothermal, ThermalConductivityModel)
+{
+private:
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+public:
+    typedef ThermalConductivitySomerton<Scalar, Indices> type;
+};
+
+//! temperature is already written by the isothermal model
+SET_BOOL_PROP(RichardsNI, NiOutputLevel, 0);
+
+//////////////////////////////////////////////////////////////////
+// Property values for isothermal model required for the general non-isothermal model
+//////////////////////////////////////////////////////////////////
+
+// set isothermal Model
+SET_TYPE_PROP(Richards, IsothermalModel, RichardsModel<TypeTag>);
+
+// set isothermal FluxVariables
+SET_TYPE_PROP(Richards, IsothermalFluxVariables, RichardsFluxVariables<TypeTag>);
+
+//set isothermal VolumeVariables
+SET_TYPE_PROP(Richards, IsothermalVolumeVariables, RichardsVolumeVariables<TypeTag>);
+
+//set isothermal LocalResidual
+SET_TYPE_PROP(Richards, IsothermalLocalResidual, RichardsLocalResidual<TypeTag>);
+
+//set isothermal Indices
+SET_TYPE_PROP(Richards, IsothermalIndices, RichardsIndices<TypeTag>);
+
+//set isothermal NumEq
+SET_INT_PROP(Richards, IsothermalNumEq, 1);
+
 // \}
 }
 
diff --git a/dumux/implicit/richards/richardsvolumevariables.hh b/dumux/implicit/richards/richardsvolumevariables.hh
index 0e24a53f45cebcefebd4bff587b580594c3fb428..0668e6e379f97ad60d03ffe4a7eec61d9c8c8873 100644
--- a/dumux/implicit/richards/richardsvolumevariables.hh
+++ b/dumux/implicit/richards/richardsvolumevariables.hh
@@ -140,6 +140,11 @@ public:
 
         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));
+
     }
 
     /*!
@@ -251,6 +256,14 @@ protected:
     {
         return problem.temperatureAtPos(fvGeometry.subContVol[scvIdx].global);
     }
+    
+    template<class ParameterCache>
+    static Scalar enthalpy_(const FluidState& fluidState,
+                            const ParameterCache& paramCache,
+                            int phaseIdx)
+    {
+        return 0;
+    }
 
     /*!
      * \brief Called by update() to compute the energy related quantities
diff --git a/test/implicit/1p/1pnispatialparams.hh b/test/implicit/1p/1pnispatialparams.hh
index 3a0540d9d38d1ae7843c9ff723f78006c1296b1d..b78135c0e1b98566cef97d5d193293e8bf83bd3f 100644
--- a/test/implicit/1p/1pnispatialparams.hh
+++ b/test/implicit/1p/1pnispatialparams.hh
@@ -144,7 +144,7 @@ public:
     }
 
     /*!
-     * \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
+     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
      *
      * \param element The finite element
      * \param fvGeometry The finite volume geometry
diff --git a/test/implicit/1p2c/1p2cnispatialparams.hh b/test/implicit/1p2c/1p2cnispatialparams.hh
index 8fe220f39884ab9062e5f5fe6ca8f81a51bf926a..2a729aae4767465621d1dd3012700527e709950d 100644
--- a/test/implicit/1p2c/1p2cnispatialparams.hh
+++ b/test/implicit/1p2c/1p2cnispatialparams.hh
@@ -144,7 +144,7 @@ public:
     }
 
     /*!
-     * \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
+     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
      *
      * \param element The finite element
      * \param fvGeometry The finite volume geometry
diff --git a/test/implicit/1p2c/1p2coutflowspatialparams.hh b/test/implicit/1p2c/1p2coutflowspatialparams.hh
index f12a3d5ea772951fc9797a42db949136165b6a08..42259673d93a3b4dd0e8324306e7163cf633f826 100644
--- a/test/implicit/1p2c/1p2coutflowspatialparams.hh
+++ b/test/implicit/1p2c/1p2coutflowspatialparams.hh
@@ -146,7 +146,7 @@ public:
     }
 
     /*!
-     * \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
+     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
      *
      * \param element The finite element
      * \param fvGeometry The finite volume geometry
diff --git a/test/implicit/2p2c/injectionspatialparams.hh b/test/implicit/2p2c/injectionspatialparams.hh
index beb335273d9ecd956d60383cf5cbbf0add126748..0ca52f7a94241cb4469e397354fd5c3fa7991aaa 100644
--- a/test/implicit/2p2c/injectionspatialparams.hh
+++ b/test/implicit/2p2c/injectionspatialparams.hh
@@ -200,7 +200,7 @@ public:
     }
 
     /*!
-     * \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the solid
+     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the solid
      *
      * This is only required for non-isothermal models.
      *
diff --git a/test/implicit/2p2c/waterairspatialparams.hh b/test/implicit/2p2c/waterairspatialparams.hh
index d7a42d1a751b6ba7dfea6a010609901dc0716fbf..95d02588ef8de8566b43be68280f90409e00bfbb 100644
--- a/test/implicit/2p2c/waterairspatialparams.hh
+++ b/test/implicit/2p2c/waterairspatialparams.hh
@@ -195,7 +195,7 @@ public:
     }
 
     /*!
-     * \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
+     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
      *
      * \param element The finite element
      * \param fvGeometry The finite volume geometry
diff --git a/test/implicit/3p/3pnispatialparams.hh b/test/implicit/3p/3pnispatialparams.hh
index 712904fedb99023bc75c33ae6ad8915728c3707a..a7fa4babf09994834b3a617a01b4df01adc11423 100644
--- a/test/implicit/3p/3pnispatialparams.hh
+++ b/test/implicit/3p/3pnispatialparams.hh
@@ -187,7 +187,7 @@ public:
     }
 
     /*!
-     * \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
+     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
      *
      * \param element The finite element
      * \param fvGeometry The finite volume geometry
diff --git a/test/implicit/3p3c/columnxylolspatialparams.hh b/test/implicit/3p3c/columnxylolspatialparams.hh
index 761b0cc591f96b1cc64dd6ab1de539dca34c3a17..3ac489d1c35fafeabbdb637642e83caba934d007 100644
--- a/test/implicit/3p3c/columnxylolspatialparams.hh
+++ b/test/implicit/3p3c/columnxylolspatialparams.hh
@@ -235,7 +235,7 @@ public:
 
 
     /*!
-     * \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
+     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
      *
      * \param element The finite element
      * \param fvGeometry The finite volume geometry
diff --git a/test/implicit/3p3c/infiltrationspatialparameters.hh b/test/implicit/3p3c/infiltrationspatialparameters.hh
index bf12030d4067830316ca2ab7d7c76b7fb4c59edb..bc6315582955803d2de374630024289887731fe0 100644
--- a/test/implicit/3p3c/infiltrationspatialparameters.hh
+++ b/test/implicit/3p3c/infiltrationspatialparameters.hh
@@ -206,55 +206,11 @@ public:
             * (1 - porosity(element, fvGeometry, scvIdx));
     }
 
-    /*!
-     * \brief Calculate the heat flux \f$[W/m^2]\f$ through the
-     *        rock matrix based on the temperature gradient \f$[K / m]\f$
-     *
-     * This is only required for non-isothermal models.
-     *
-     * \param heatFlux The resulting heat flux vector
-     * \param fluxVars The flux variables
-     * \param elemVolVars The volume variables
-     * \param tempGrad The temperature gradient
-     * \param element The current finite element
-     * \param fvGeometry The finite volume geometry of the current element
-     * \param scvfIdx The local index of the sub-control volume face where
-     *                    the matrix heat flux should be calculated
-     */
-    void matrixHeatFlux(DimVector &heatFlux,
-                        const FluxVariables &fluxVars,
-                        const ElementVolumeVariables &elemVolVars,
-                        const DimVector &tempGrad,
-                        const Element &element,
-                        const FVElementGeometry &fvGeometry,
-                        int scvfIdx) const
-    {
-        static const Scalar lDry = 0.35;
-        static const Scalar lsw1 = 1.8;
-        static const Scalar lsn1 = 0.65;
-
-        // arithmetic mean of the liquid saturation and the porosity
-        const int i = fluxVars.face().i;
-        const int j = fluxVars.face().j;
-        Scalar sw = std::max(0.0, (elemVolVars[i].saturation(wPhaseIdx) +
-                                   elemVolVars[j].saturation(wPhaseIdx)) / 2);
-        Scalar sn = std::max(0.0, (elemVolVars[i].saturation(nPhaseIdx) +
-                                   elemVolVars[j].saturation(nPhaseIdx)) / 2);
-
-        // the heat conductivity of the matrix. in general this is a
-        // tensorial value, but we assume isotropic heat conductivity.
-        Scalar heatCond = lDry + sqrt(sw) * (lsw1-lDry) + sqrt(sn) * (lsn1-lDry);
-
-        // the matrix heat flux is the negative temperature gradient
-        // times the heat conductivity.
-        heatFlux = tempGrad;
-        heatFlux *= -heatCond;
-    }
-
     const MaterialLawParams& materialLawParams() const
     {
         return MaterialParams_;
     }
+
 private:
     bool isFineMaterial_(const GlobalPosition &globalPos) const
     { return
diff --git a/test/implicit/3p3c/kuevettespatialparams.hh b/test/implicit/3p3c/kuevettespatialparams.hh
index 47689070df69a7e576beb2f8e6d027ec971608f2..14f8cccddc9164b5dcc34dd09d9d6924afebdfe1 100644
--- a/test/implicit/3p3c/kuevettespatialparams.hh
+++ b/test/implicit/3p3c/kuevettespatialparams.hh
@@ -226,52 +226,7 @@ public:
     }
 
     /*!
-     * \brief Calculate the heat flux \f$[W/m^2]\f$ through the
-     *        rock matrix based on the temperature gradient \f$[K / m]\f$
-     *
-     * This is only required for non-isothermal models.
-     *
-     * \param heatFlux The resulting heat flux vector
-     * \param fluxVars The flux variables
-     * \param elemVolVars The volume variables
-     * \param tempGrad The temperature gradient
-     * \param element The current finite element
-     * \param fvGeometry The finite volume geometry of the current element
-     * \param fIdx The local index of the sub-control volume face where
-     *                    the matrix heat flux should be calculated
-     */
-    void matrixHeatFlux(DimVector &heatFlux,
-                        const FluxVariables &fluxVars,
-                        const ElementVolumeVariables &elemVolVars,
-                        const DimVector &tempGrad,
-                        const Element &element,
-                        const FVElementGeometry &fvGeometry,
-                        const int fIdx) const
-    {
-        static const Scalar lDry = 0.35;
-        static const Scalar lsw1 = 1.8;
-        static const Scalar lsn1 = 0.65;
-
-        // arithmetic mean of the liquid saturation and the porosity
-        const int i = fluxVars.face().i;
-        const int j = fluxVars.face().j;
-        Scalar sw = std::max(0.0, (elemVolVars[i].saturation(wPhaseIdx) +
-                                   elemVolVars[j].saturation(wPhaseIdx)) / 2);
-        Scalar sn = std::max(0.0, (elemVolVars[i].saturation(nPhaseIdx) +
-                                   elemVolVars[j].saturation(nPhaseIdx)) / 2);
-
-        // the heat conductivity of the matrix. in general this is a
-        // tensorial value, but we assume isotropic heat conductivity.
-        Scalar heatCond = lDry + sqrt(sw) * (lsw1-lDry) + sqrt(sn) * (lsn1-lDry);
-
-        // the matrix heat flux is the negative temperature gradient
-        // times the heat conductivity.
-        heatFlux = tempGrad;
-        heatFlux *= -heatCond;
-    }
-
-    /*!
-     * \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the porous material.
+     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
      *
      * \param element The finite element
      * \param fvGeometry The finite volume geometry
diff --git a/test/implicit/co2/heterogeneousspatialparameters.hh b/test/implicit/co2/heterogeneousspatialparameters.hh
index 92c466826226a34dc79c8e4fe7cd24895c745232..d0480c8b94893016c93742b5b5037b03063ab117 100644
--- a/test/implicit/co2/heterogeneousspatialparameters.hh
+++ b/test/implicit/co2/heterogeneousspatialparameters.hh
@@ -231,7 +231,7 @@ public:
     }
 
     /*!
-     * \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the solid
+     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the solid
      *
      * This is only required for non-isothermal models.
      *
diff --git a/test/implicit/richards/CMakeLists.txt b/test/implicit/richards/CMakeLists.txt
index 8defc79b558fb479de3c9f0818a2201aa8f56c23..354b7ad94b1ce96d1b8de80286a48aa820bb357c 100644
--- a/test/implicit/richards/CMakeLists.txt
+++ b/test/implicit/richards/CMakeLists.txt
@@ -6,6 +6,7 @@ if(MPI_CXX_FOUND)
     ${CMAKE_CURRENT_BINARY_DIR}/s0002-p0000-richardslensbox-00008.vtu
     "${MPIEXEC} -np 2 ${CMAKE_CURRENT_BINARY_DIR}/test_boxrichards")
 else()
+# isothermal tests
   add_dumux_test(test_boxrichards test_boxrichards test_boxrichards.cc
     ${CMAKE_SOURCE_DIR}/bin/runTest.sh
     ${CMAKE_SOURCE_DIR}/bin/fuzzycomparevtu.py
@@ -20,3 +21,32 @@ add_dumux_test(test_ccrichards test_ccrichards test_ccrichards.cc
   ${CMAKE_SOURCE_DIR}/test/references/richardslenscc-reference.vtu
   ${CMAKE_CURRENT_BINARY_DIR}/richardslenscc-00008.vtu
   ${CMAKE_CURRENT_BINARY_DIR}/test_ccrichards)
+
+# non-isothermal tests
+add_dumux_test(test_boxrichardsniconvection test_boxrichardsniconvection test_boxrichardsniconvection.cc
+  ${CMAKE_SOURCE_DIR}/bin/runTest.sh
+  ${CMAKE_SOURCE_DIR}/bin/fuzzycomparevtu.py
+  ${CMAKE_SOURCE_DIR}/test/references/richardsniconvectionbox-reference.vtu
+  ${CMAKE_CURRENT_BINARY_DIR}/test_boxrichardsniconvection-00011.vtu
+  ${CMAKE_CURRENT_BINARY_DIR}/test_boxrichardsniconvection)
+
+add_dumux_test(test_ccrichardsniconvection test_ccrichardsniconvection test_ccrichardsniconvection.cc
+  ${CMAKE_SOURCE_DIR}/bin/runTest.sh
+  ${CMAKE_SOURCE_DIR}/bin/fuzzycomparevtu.py
+  ${CMAKE_SOURCE_DIR}/test/references/richardsniconvectioncc-reference.vtu
+  ${CMAKE_CURRENT_BINARY_DIR}/test_ccrichardsniconvection-00011.vtu
+  ${CMAKE_CURRENT_BINARY_DIR}/test_ccrichardsniconvection)
+
+add_dumux_test(test_boxrichardsniconduction test_boxrichardsniconduction test_boxrichardsniconduction.cc
+  ${CMAKE_SOURCE_DIR}/bin/runTest.sh
+  ${CMAKE_SOURCE_DIR}/bin/fuzzycomparevtu.py
+  ${CMAKE_SOURCE_DIR}/test/references/richardsniconductionbox-reference.vtu
+  ${CMAKE_CURRENT_BINARY_DIR}/test_boxrichardsniconduction-00007.vtu
+  ${CMAKE_CURRENT_BINARY_DIR}/test_boxrichardsniconduction)
+
+add_dumux_test(test_ccrichardsniconduction test_ccrichardsniconduction test_ccrichardsniconduction.cc
+  ${CMAKE_SOURCE_DIR}/bin/runTest.sh
+  ${CMAKE_SOURCE_DIR}/bin/fuzzycomparevtu.py
+  ${CMAKE_SOURCE_DIR}/test/references/richardsniconductioncc-reference.vtu
+  ${CMAKE_CURRENT_BINARY_DIR}/test_ccrichardsniconduction-00007.vtu
+  ${CMAKE_CURRENT_BINARY_DIR}/test_ccrichardsniconduction)
\ No newline at end of file
diff --git a/test/implicit/richards/Makefile.am b/test/implicit/richards/Makefile.am
index 20e21e07a2c83a7db4a1ad7a73f4f9a57f7b7cd6..4bb2a4e7d759aa159befa9f16eefce9a85e9226b 100644
--- a/test/implicit/richards/Makefile.am
+++ b/test/implicit/richards/Makefile.am
@@ -1,9 +1,19 @@
-check_PROGRAMS = test_boxrichards test_ccrichards
+check_PROGRAMS = \
+	test_boxrichards \
+	test_ccrichards \
+	test_boxrichardsniconvection \ 
+	test_ccrichardsniconvection \ 
+	test_boxrichardsniconduction \
+	test_ccrichardsniconduction
 
 noinst_HEADERS := $(wildcard *.hh)
 EXTRA_DIST:=$(wildcard *.input) $(wildcard grids/*.dgf) CMakeLists.txt
 
 test_boxrichards_SOURCES = test_boxrichards.cc 
 test_ccrichards_SOURCES = test_ccrichards.cc 
+test_boxrichardsniconvection_SOURCES = test_boxrichardsniconvection.cc
+test_ccrichardsniconvection_SOURCES = test_ccrichardsniconvection.cc
+test_boxrichardsniconduction_SOURCES = test_boxrichardsniconduction.cc
+test_ccrichardsniconduction_SOURCES = test_ccrichardsniconduction.cc
 
 include $(top_srcdir)/am/global-rules
diff --git a/test/implicit/richards/richardsniconductionproblem.hh b/test/implicit/richards/richardsniconductionproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..b4b6039859dca74bc1e34b907fc733396726e12f
--- /dev/null
+++ b/test/implicit/richards/richardsniconductionproblem.hh
@@ -0,0 +1,363 @@
+// -*- 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 Test for the RichardsModel in combination with the NI model for a conduction problem:
+ * The simulation domain is a tube where with an elevated temperature on the left hand side.
+ */
+#ifndef DUMUX_RICHARDS_CONDUCTION_PROBLEM_HH
+#define DUMUX_RICHARDS_CONDUCTION_PROBLEM_HH
+
+#include <math.h>
+#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
+
+#include <dumux/implicit/common/implicitporousmediaproblem.hh>
+#include <dumux/implicit/richards/richardsmodel.hh>
+#include <dumux/implicit/richards/richardsproblem.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
+#include <dumux/material/fluidsystems/h2on2fluidsystem.hh>
+#include "richardsnispatialparams.hh"
+
+namespace Dumux
+{
+
+template <class TypeTag>
+class RichardsNIConductionProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(RichardsNIConductionProblem, INHERITS_FROM(RichardsNI, RichardsNISpatialParams));
+NEW_TYPE_TAG(RichardsNIConductionBoxProblem, INHERITS_FROM(BoxModel, RichardsNIConductionProblem));
+NEW_TYPE_TAG(RichardsNIConductionCCProblem, INHERITS_FROM(CCModel, RichardsNIConductionProblem));
+
+// Set the grid type
+SET_TYPE_PROP(RichardsNIConductionProblem, Grid, Dune::YaspGrid<2>);
+
+// Set the problem property
+SET_TYPE_PROP(RichardsNIConductionProblem, Problem,
+              Dumux::RichardsNIConductionProblem<TypeTag>);
+
+// Set the fluid system
+SET_TYPE_PROP(RichardsNIConductionProblem, FluidSystem, FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), false>);
+
+// Set the spatial parameters
+SET_TYPE_PROP(RichardsNIConductionProblem,
+              SpatialParams,
+              Dumux::RichardsNISpatialParams<TypeTag>);
+
+// Enable velocity output
+SET_BOOL_PROP(RichardsNIConductionProblem, VtkAddVelocity, true);
+
+// Disable gravity
+SET_BOOL_PROP(RichardsNIConductionProblem, ProblemEnableGravity, false);
+}
+
+/*!
+ * \ingroup RichardsModel
+ * \ingroup ImplicitTestProblems
+ *
+ * \brief Test for the RichardsModel in combination with the NI model for a conduction problem:
+ * The simulation domain is a tube where with an elevated temperature on the left hand side.
+ *
+ * Initially the domain is fully saturated with water at a constant temperature.
+ * On the left hand side there is a Dirichlet boundary condition with an increased temperature and on the right hand side
+ * a Dirichlet boundary with constant pressure, saturation and temperature is applied.
+ *
+ * The results are compared to an analytical solution for a diffusion process:
+  \f[
+     T =T_{high} + (T_{init} - T_{high})erf \left(0.5\sqrt{\frac{x^2 S_{total}}{t \lambda_{eff}}}\right)
+ \f]
+ *
+ * This problem uses the \ref RichardsModel and \ref NIModel model.
+ *
+ * To run the simulation execute the following line in shell: <br>
+ * <tt>./test_boxrichardsniconduction -ParameterFile ./test_boxrichardsniconduction.input</tt> or <br>
+ * <tt>./test_ccrichardsniconduction -ParameterFile ./test_ccrichardsniconduction.input</tt>
+ */
+template <class TypeTag>
+class RichardsNIConductionProblem : public RichardsBoxProblem<TypeTag>
+{
+    typedef RichardsBoxProblem<TypeTag> ParentType;
+
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+    typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+    typedef typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel) ThermalConductivityModel;
+    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+
+    // copy some indices for convenience
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    enum {
+        // world dimension
+        dimWorld = GridView::dimensionworld
+    };
+
+    enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+    enum { dofCodim = isBox ? dimWorld : 0 };
+
+    enum {
+        // indices of the primary variables
+        pressureIdx = Indices::pwIdx,
+        wPhaseIdx = Indices::wPhaseIdx,
+        nPhaseIdx = Indices::nPhaseIdx,
+        temperatureIdx = Indices::temperatureIdx
+    };
+    enum {
+        // index of the transport equation
+        contiEqIdx = Indices::contiEqIdx,
+        energyEqIdx = Indices::energyEqIdx
+    };
+
+    typedef typename GridView::template Codim<0>::Entity Element;
+    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+    typedef typename GridView::Intersection Intersection;
+
+    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+    RichardsNIConductionProblem(TimeManager &timeManager, const GridView &gridView)
+        : ParentType(timeManager, gridView)
+        , eps_(1e-6)
+    {
+        //initialize fluid system
+        FluidSystem::init();
+
+        name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
+        outputInterval_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Problem, OutputInterval);
+        temperatureHigh_ = 300.;
+    }
+
+
+    bool shouldWriteOutput() const
+    {
+        return
+            this->timeManager().timeStepIndex() == 0 ||
+            this->timeManager().timeStepIndex() % outputInterval_ == 0 ||
+            this->timeManager().episodeWillBeOver() ||
+            this->timeManager().willBeFinished();
+    }
+
+    /*!
+     * \brief Append all quantities of interest which can be derived
+     *        from the solution of the current time step to the VTK
+     *        writer.
+     */
+    void addOutputVtkFields()
+    {
+        //Here we calculate the analytical solution
+        typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
+        unsigned numDofs = this->model().numDofs();
+
+        //create required scalar fields
+        ScalarField *temperatureExact = this->resultWriter().allocateManagedBuffer(numDofs);
+
+        FVElementGeometry fvGeometry;
+        VolumeVariables volVars;
+
+        ElementIterator eIt = this->gridView().template begin<0>();
+        ElementIterator eEndIt = this->gridView().template end<0>();
+        fvGeometry.update(this->gridView(), *eIt);
+        PrimaryVariables initialPriVars(0);
+        GlobalPosition globalPos(0);
+        initial_(initialPriVars, globalPos);
+
+        //update the constant volume variables
+        volVars.update(initialPriVars, *this, *eIt, fvGeometry, 0, false);
+
+        Scalar porosity = this->spatialParams().porosity(*eIt, fvGeometry, 0);
+        Scalar densityW = volVars.density(wPhaseIdx);
+        Scalar heatCapacityW = FluidSystem::heatCapacity(volVars.fluidState(), 0);
+        Scalar effectiveHeatCapacityS = this->spatialParams().heatCapacity(*eIt, fvGeometry, 0);
+        Scalar storage = densityW*heatCapacityW*porosity + effectiveHeatCapacityS;
+        Scalar effectiveThermalConductivity = ThermalConductivityModel::effectiveThermalConductivity(volVars, this->spatialParams(),
+                *eIt, fvGeometry, 0);
+        Scalar time = std::max(this->timeManager().time() + this->timeManager().timeStepSize(), 1e-10);
+
+        for (; eIt != eEndIt; ++eIt)
+        {
+            fvGeometry.update(this->gridView(), *eIt);
+            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+            {
+                int globalIdx = this->model().dofMapper().map(*eIt, scvIdx, dofCodim);
+                if (isBox)
+                    globalPos = eIt->geometry().corner(scvIdx);
+                else
+                    globalPos = eIt->geometry().center();
+
+                (*temperatureExact)[globalIdx] = temperatureHigh_ + (initialPriVars[temperatureIdx] - temperatureHigh_)
+                                                             *std::erf(0.5*std::sqrt(globalPos[0]*globalPos[0]*storage/time/effectiveThermalConductivity));
+            }
+        }
+        this->resultWriter().attachDofData(*temperatureExact, "temperatureExact", isBox);
+    }
+    /*!
+     * \name Problem parameters
+     */
+    // \{
+
+    /*!
+     * \brief The problem name.
+     *
+     * This is used as a prefix for files generated by the simulation.
+     */
+    const char *name() const
+    {
+        return name_.c_str();
+    }
+    // \}
+
+    /*!
+     * \name Boundary conditions
+     */
+    // \{
+
+    /*!
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary segment.
+     *
+     * \param values The boundary types for the conservation equations
+     * \param globalPos The position for which the bc type should be evaluated
+     */
+    void boundaryTypesAtPos(BoundaryTypes &values, 
+                            const GlobalPosition &globalPos) const
+    {
+        if(globalPos[0] < eps_ || globalPos[0] > this->bBoxMax()[0] - eps_)
+        {
+            values.setAllDirichlet();
+        }
+        else
+        {
+            values.setAllNeumann();
+        }
+    }
+
+    /*!
+     * \brief Evaluate the boundary conditions for a dirichlet
+     *        boundary segment.
+     *
+     * \param values The dirichlet values for the primary variables
+     * \param globalPos The position for which the bc type should be evaluated
+     *
+     * For this method, the \a values parameter stores primary variables.
+     */
+    void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+    {
+        initial_(values, globalPos);
+
+        if (globalPos[0] < eps_)
+        {
+            values[temperatureIdx] = temperatureHigh_;
+        }
+    }
+
+    /*!
+     * \brief Evaluate the boundary conditions for a Neumann
+     *        boundary segment.
+     *
+     * For this method, the \a priVars parameter stores the mass flux
+     * in normal direction of each component. Negative values mean
+     * influx.
+     */
+    void neumann(PrimaryVariables &priVars,
+                 const Element &element,
+                 const FVElementGeometry &fvGeometry,
+                 const Intersection &intersection,
+                 const int scvIdx,
+                 const int boundaryFaceIdx) const
+    {
+        priVars = 0;
+    }
+
+    // \}
+
+    /*!
+     * \name Volume terms
+     */
+    // \{
+
+    /*!
+     * \brief Evaluate the source term for all phases within a given
+     *        sub-control-volume.
+     *
+     * For this method, the \a priVars parameter stores the rate mass
+     * of a component is generated or annihilate per volume
+     * unit. Positive values mean that mass is created, negative ones
+     * mean that it vanishes.
+     */
+    void sourceAtPos(PrimaryVariables &priVars,
+                     const GlobalPosition &globalPos) const
+    {
+        priVars = Scalar(0.0);
+    }
+
+    /*!
+     * \brief Returns the reference pressure [Pa] of the non-wetting
+     *        fluid phase within a finite volume
+     *
+     * This problem assumes a constant reference pressure of 1 bar.
+     *
+     * \param element The DUNE Codim<0> entity which intersects with
+     *                the finite volume in question
+     * \param fvGeometry The finite volume geometry of the element
+     * \param scvIdx The sub control volume index inside the finite
+     *               volume geometry
+     */
+    Scalar referencePressure(const Element &element,
+                             const FVElementGeometry &fvGeometry,
+                             const int scvIdx) const
+    { return 1e5; };
+
+    /*!
+     * \brief Evaluate the initial value for a control volume.
+     *
+     * \param values The initial values for the primary variables
+     * \param globalPos The position for which the initial condition should be evaluated
+     *
+     * For this method, the \a values parameter stores primary
+     * variables.
+     */
+    void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+    {
+        initial_(values, globalPos);
+    }
+
+    // \}
+
+private:
+    // the internal method for the initial condition
+    void initial_(PrimaryVariables &priVars,
+                  const GlobalPosition &globalPos) const
+    {
+        priVars[pressureIdx] = 1e5; // initial condition for the pressure
+        priVars[temperatureIdx] = 290.;
+    }
+
+    Scalar temperatureHigh_;
+    const Scalar eps_;
+    std::string name_;
+    int outputInterval_;
+};
+
+} //end namespace
+#endif // DUMUX_RICHARDSNINI_CONDUCTION_PROBLEM_HH
diff --git a/test/implicit/richards/richardsniconvectionproblem.hh b/test/implicit/richards/richardsniconvectionproblem.hh
new file mode 100644
index 0000000000000000000000000000000000000000..7d6e21a26d6020ff3a183d2fa808d9ac104220a0
--- /dev/null
+++ b/test/implicit/richards/richardsniconvectionproblem.hh
@@ -0,0 +1,396 @@
+// -*- 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 Test for the RichardsModel in combination with the NI model for a convection problem:
+ * The simulation domain is a tube where water with an elevated temperature is injected
+ * at a constant rate on the left hand side.
+ */
+#ifndef DUMUX_1PNI_CONVECTION_PROBLEM_HH
+#define DUMUX_1PNI_CONVECTION_PROBLEM_HH
+
+#include <math.h>
+#include <dune/grid/io/file/dgfparser/dgfyasp.hh>
+
+#include <dumux/implicit/common/implicitporousmediaproblem.hh>
+#include <dumux/implicit/richards/richardsmodel.hh>
+#include <dumux/implicit/richards/richardsproblem.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh>
+#include <dumux/material/fluidsystems/h2on2fluidsystem.hh>
+#include "richardsnispatialparams.hh"
+
+namespace Dumux
+{
+
+template <class TypeTag>
+class RichardsNIConvectionProblem;
+
+namespace Properties
+{
+NEW_TYPE_TAG(RichardsNIConvectionProblem, INHERITS_FROM(RichardsNI, RichardsNISpatialParams));
+NEW_TYPE_TAG(RichardsNIConvectionBoxProblem, INHERITS_FROM(BoxModel, RichardsNIConvectionProblem));
+NEW_TYPE_TAG(RichardsNIConvectionCCProblem, INHERITS_FROM(CCModel, RichardsNIConvectionProblem));
+
+// Set the grid type
+SET_TYPE_PROP(RichardsNIConvectionProblem, Grid, Dune::YaspGrid<2>);
+
+// Set the problem property
+SET_TYPE_PROP(RichardsNIConvectionProblem, Problem,
+              Dumux::RichardsNIConvectionProblem<TypeTag>);
+
+// Set the fluid system
+SET_TYPE_PROP(RichardsNIConvectionProblem, FluidSystem, FluidSystems::H2ON2<typename GET_PROP_TYPE(TypeTag, Scalar), false>);
+
+// Set the spatial parameters
+SET_TYPE_PROP(RichardsNIConvectionProblem,
+              SpatialParams,
+              Dumux::RichardsNISpatialParams<TypeTag>);
+
+// Enable velocity output
+SET_BOOL_PROP(RichardsNIConvectionProblem, VtkAddVelocity, true);
+
+// Disable gravity
+SET_BOOL_PROP(RichardsNIConvectionProblem, ProblemEnableGravity, false);
+}
+
+
+/*!
+ * \ingroup RichardsModel
+ * \ingroup ImplicitTestProblems
+ *
+ * \brief Test for the RichardsModel in combination with the NI model for a convection problem:
+ * The simulation domain is a tube where water with an elevated temperature is injected
+ * at a constant rate on the left hand side.
+ * 
+ * Initially the domain is fully saturated with water at a constant temperature.
+ * On the left hand side water is injected at a constant rate and on the right hand side
+ * a Dirichlet boundary with constant pressure, saturation and temperature is applied.
+ *
+ * The results are compared to an analytical solution where a retarded front velocity is calculated as follows:
+  \f[
+     v_{Front}=\frac{q S_{water}}{\phi S_{total}}
+ \f]
+ *
+ * The result of the analytical solution is written into the vtu files.
+ * This problem uses the \ref RichardsModel and \ref NIModel model.
+ *
+ * To run the simulation execute the following line in shell: <br>
+ * <tt>./test_boxrichardsniconvection -ParameterFile ./test_boxrichardsniconvection.input</tt> or <br>
+ * <tt>./test_ccrichardsniconvection -ParameterFile ./test_ccrichardsniconvection.input</tt>
+ */
+template <class TypeTag>
+class RichardsNIConvectionProblem : public RichardsBoxProblem<TypeTag>
+{
+    typedef RichardsBoxProblem<TypeTag> ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+    typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
+    typedef typename GET_PROP_TYPE(TypeTag, BoundaryTypes) BoundaryTypes;
+    typedef typename GET_PROP_TYPE(TypeTag, TimeManager) TimeManager;
+    typedef typename GET_PROP_TYPE(TypeTag, ThermalConductivityModel) ThermalConductivityModel;
+    typedef typename GET_PROP_TYPE(TypeTag, ElementVolumeVariables) ElementVolumeVariables;
+    typedef typename GET_PROP_TYPE(TypeTag, VolumeVariables) VolumeVariables;
+    typedef Dumux::H2O<Scalar> IapwsH2O;
+
+    // copy some indices for convenience
+    typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
+    enum {
+        // world dimension
+        dimWorld = GridView::dimensionworld
+    };
+
+    enum { isBox = GET_PROP_VALUE(TypeTag, ImplicitIsBox) };
+    enum { dofCodim = isBox ? dimWorld : 0 };
+
+    enum {
+        // indices of the primary variables
+        pressureIdx = Indices::pwIdx,
+        wPhaseIdx = Indices::wPhaseIdx,
+        nPhaseIdx = Indices::nPhaseIdx,
+        temperatureIdx = Indices::temperatureIdx
+    };
+    enum {
+        // index of the transport equation
+        contiEqIdx = Indices::contiEqIdx,
+        energyEqIdx = Indices::energyEqIdx
+    };
+
+    typedef typename GridView::template Codim<0>::Entity Element;
+    typedef typename GridView::template Codim<0>::Iterator ElementIterator;
+    typedef typename GridView::Intersection Intersection;
+
+    typedef Dune::FieldVector<Scalar, dimWorld> GlobalPosition;
+
+public:
+    RichardsNIConvectionProblem(TimeManager &timeManager, const GridView &gridView)
+        : ParentType(timeManager, gridView)
+        , eps_(1e-6)
+    {
+        //initialize fluid system
+        FluidSystem::init();
+
+        name_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, std::string, Problem, Name);
+        outputInterval_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, int, Problem, OutputInterval);
+        darcyVelocity_ = GET_RUNTIME_PARAM_FROM_GROUP(TypeTag, Scalar, Problem, DarcyVelocity);
+        temperatureHigh_ = 291.;
+        temperatureLow_ = 290.;
+        pressureHigh_ = 2e5;
+        pressureLow_ = 1e5;
+    }
+
+
+    bool shouldWriteOutput() const
+    {
+        return
+            this->timeManager().timeStepIndex() == 0 ||
+            this->timeManager().timeStepIndex() % outputInterval_ == 0 ||
+            this->timeManager().episodeWillBeOver() ||
+            this->timeManager().willBeFinished();
+    }
+
+    /*!
+     * \brief Append all quantities of interest which can be derived
+     *        from the solution of the current time step to the VTK
+     *        writer.
+     */
+    void addOutputVtkFields()
+    {
+        //Here we calculate the analytical solution
+        typedef Dune::BlockVector<Dune::FieldVector<double, 1> > ScalarField;
+        unsigned numDofs = this->model().numDofs();
+
+        //create required scalar fields
+        ScalarField *temperatureExact = this->resultWriter().allocateManagedBuffer(numDofs);
+
+        FVElementGeometry fvGeometry;
+        VolumeVariables volVars;
+
+        ElementIterator eIt = this->gridView().template begin<0>();
+        ElementIterator eEndIt = this->gridView().template end<0>();
+        fvGeometry.update(this->gridView(), *eIt);
+        PrimaryVariables initialPriVars(0);
+        GlobalPosition globalPos(0);
+        initial_(initialPriVars, globalPos);
+
+        //update the constant volume variables
+        volVars.update(initialPriVars, *this, *eIt, fvGeometry, 0, false);
+
+        Scalar porosity = this->spatialParams().porosity(*eIt, fvGeometry, 0);
+        Scalar densityW = volVars.density(wPhaseIdx);
+        Scalar heatCapacityW = FluidSystem::heatCapacity(volVars.fluidState(), 0);
+        Scalar effectiveHeatCapacityS = this->spatialParams().heatCapacity(*eIt, fvGeometry, 0);
+        Scalar storageW =  densityW*heatCapacityW*porosity;
+        Scalar storageTotal = storageW + effectiveHeatCapacityS;
+        std::cout<<"storage: "<<storageTotal<<std::endl;
+        Scalar time = std::max(this->timeManager().time() + this->timeManager().timeStepSize(), 1e-10);
+        Scalar retardedFrontVelocity = darcyVelocity_*storageW/storageTotal/porosity;
+        std::cout<<"retarded velocity: "<<retardedFrontVelocity<<std::endl;
+
+        for (; eIt != eEndIt; ++eIt)
+        {
+            fvGeometry.update(this->gridView(), *eIt);
+            for (int scvIdx = 0; scvIdx < fvGeometry.numScv; ++scvIdx)
+            {
+                int globalIdx = this->model().dofMapper().map(*eIt, scvIdx, dofCodim);
+                if (isBox)
+                    globalPos = eIt->geometry().corner(scvIdx);
+                else
+                    globalPos = eIt->geometry().center();
+
+                (*temperatureExact)[globalIdx] = globalPos[0] < retardedFrontVelocity*time ? temperatureHigh_ : temperatureLow_;
+            }
+        }
+        this->resultWriter().attachDofData(*temperatureExact, "temperatureExact", isBox);
+    }
+    /*!
+     * \name Problem parameters
+     */
+    // \{
+
+    /*!
+     * \brief The problem name.
+     *
+     * This is used as a prefix for files generated by the simulation.
+     */
+    const char *name() const
+    {
+        return name_.c_str();
+    }
+
+    // \}
+
+    /*!
+     * \name Boundary conditions
+     */
+    // \{
+
+    /*!
+     * \brief Specifies which kind of boundary condition should be
+     *        used for which equation on a given boundary segment.
+     *
+     * \param values The boundary types for the conservation equations
+     * \param globalPos The position for which the bc type should be evaluated
+     */
+    void boundaryTypesAtPos(BoundaryTypes &values, 
+                            const GlobalPosition &globalPos) const
+    {
+        if(globalPos[0] > this->bBoxMax()[0] - eps_)
+        {
+            values.setAllDirichlet();
+        }
+        else
+        {
+            values.setAllNeumann();
+        }
+    }
+
+    /*!
+     * \brief Evaluate the boundary conditions for a dirichlet
+     *        boundary segment.
+     *
+     * \param values The dirichlet values for the primary variables
+     * \param globalPos The position for which the bc type should be evaluated
+     *
+     * For this method, the \a values parameter stores primary variables.
+     */
+    void dirichletAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+    {
+        initial_(values, globalPos);
+    }
+
+    /*!
+     * \brief Evaluates the boundary conditions for a Neumann
+     *        boundary segment in dependency on the current solution.
+     *
+     * \param values Stores the Neumann values for the conservation equations in
+     *               \f$ [ \textnormal{unit of conserved quantity} / (m^(dim-1) \cdot s )] \f$
+     * \param element The finite element
+     * \param fvGeometry The finite volume geometry of the element
+     * \param intersection The intersection between element and boundary
+     * \param scvIdx The local index of the sub-control volume
+     * \param boundaryFaceIdx The index of the boundary face
+     * \param elemVolVars All volume variables for the element
+     *
+     * This method is used for cases, when the Neumann condition depends on the
+     * solution and requires some quantities that are specific to the fully-implicit method.
+     * The \a values store the mass flux of each phase normal to the boundary.
+     * Negative values indicate an inflow.
+     */
+     void solDependentNeumann(PrimaryVariables &values,
+                      const Element &element,
+                      const FVElementGeometry &fvGeometry,
+                      const Intersection &intersection,
+                      const int scvIdx,
+                      const int boundaryFaceIdx,
+                      const ElementVolumeVariables &elemVolVars) const
+    {
+        values = 0;
+        GlobalPosition globalPos(0);
+        if (isBox)
+            globalPos = fvGeometry.boundaryFace[boundaryFaceIdx].ipGlobal;
+        else
+            globalPos = fvGeometry.boundaryFace[boundaryFaceIdx].ipGlobal;
+
+        if(globalPos[0] < eps_)
+        {
+            values[pressureIdx] = -darcyVelocity_*elemVolVars[scvIdx].density(wPhaseIdx);
+            values[temperatureIdx] = -darcyVelocity_*elemVolVars[scvIdx].density(wPhaseIdx)
+                                     *IapwsH2O::liquidEnthalpy(temperatureHigh_, elemVolVars[scvIdx].pressure(wPhaseIdx));
+        }
+    }
+
+    // \}
+
+    /*!
+     * \name Volume terms
+     */
+    // \{
+
+    /*!
+     * \brief Evaluate the source term for all phases within a given
+     *        sub-control-volume.
+     *
+     * For this method, the \a priVars parameter stores the rate mass
+     * of a component is generated or annihilate per volume
+     * unit. Positive values mean that mass is created, negative ones
+     * mean that it vanishes.
+     */
+    void sourceAtPos(PrimaryVariables &priVars,
+                     const GlobalPosition &globalPos) const
+    {
+        priVars = Scalar(0.0);
+    }
+
+    /*!
+     * \brief Returns the reference pressure [Pa] of the non-wetting
+     *        fluid phase within a finite volume
+     *
+     * This problem assumes a constant reference pressure of 1 bar.
+     *
+     * \param element The DUNE Codim<0> entity which intersects with
+     *                the finite volume in question
+     * \param fvGeometry The finite volume geometry of the element
+     * \param scvIdx The sub control volume index inside the finite
+     *               volume geometry
+     */
+    Scalar referencePressure(const Element &element,
+                             const FVElementGeometry &fvGeometry,
+                             const int scvIdx) const
+    { return 1e5; };
+
+    /*!
+     * \brief Evaluate the initial value for a control volume.
+     *
+     * \param values The initial values for the primary variables
+     * \param globalPos The position for which the initial condition should be evaluated
+     *
+     * For this method, the \a values parameter stores primary
+     * variables.
+     */
+    void initialAtPos(PrimaryVariables &values, const GlobalPosition &globalPos) const
+    {
+        initial_(values, globalPos);
+    }
+
+    // \}
+
+private:
+    // the internal method for the initial condition
+    void initial_(PrimaryVariables &priVars,
+                  const GlobalPosition &globalPos) const
+    {
+        priVars[pressureIdx] = pressureLow_; // initial condition for the pressure
+        priVars[temperatureIdx] = temperatureLow_;
+    }
+
+    Scalar temperatureHigh_;
+    Scalar temperatureLow_;
+    Scalar pressureHigh_;
+    Scalar pressureLow_;
+    Scalar darcyVelocity_;
+    const Scalar eps_;
+    std::string name_;
+    int outputInterval_;
+};
+
+} //end namespace
+#endif // DUMUX_RICHARDSNI_CONVECTION_PROBLEM_HH
diff --git a/test/implicit/richards/richardsnispatialparams.hh b/test/implicit/richards/richardsnispatialparams.hh
new file mode 100644
index 0000000000000000000000000000000000000000..9e92fd47ee7da10698907e97341d3335097e3613
--- /dev/null
+++ b/test/implicit/richards/richardsnispatialparams.hh
@@ -0,0 +1,222 @@
+// -*- 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 Definition of the spatial parameters for the non-isothermal Richards problems.
+ */
+#ifndef DUMUX_RICHARDSNI_SPATIAL_PARAMS_HH
+#define DUMUX_RICHARDSNI_SPATIAL_PARAMS_HH
+
+#include <dumux/implicit/richards/richardsmodel.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/efftoabslaw.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/linearmaterial.hh>
+#include <dumux/material/fluidmatrixinteractions/2p/regularizedvangenuchten.hh>
+#include <dumux/material/spatialparams/implicitspatialparams.hh>
+
+namespace Dumux
+{
+
+/*!
+ * \ingroup RichardsModel
+ * \ingroup ImplicitTestProblems
+ *
+ * \brief Definition of the spatial parameters for the RichardsNI problems.
+ */
+
+//forward declaration
+template<class TypeTag>
+class RichardsNISpatialParams;
+
+namespace Properties
+{
+// The spatial parameters TypeTag
+NEW_TYPE_TAG(RichardsNISpatialParams);
+
+// Set the spatial parameters
+SET_TYPE_PROP(RichardsNISpatialParams, SpatialParams, Dumux::RichardsNISpatialParams<TypeTag>);
+
+// Set the material law
+SET_PROP(RichardsNISpatialParams, MaterialLaw)
+{
+private:
+    // define the material law which is parameterized by effective
+    // saturations
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+    typedef RegularizedVanGenuchten<Scalar> EffectiveLaw;
+public:
+    // define the material law parameterized by absolute saturations
+    typedef EffToAbsLaw<EffectiveLaw> type;
+};
+}
+
+
+template<class TypeTag>
+class RichardsNISpatialParams : public ImplicitSpatialParams<TypeTag>
+{
+    typedef ImplicitSpatialParams<TypeTag> ParentType;
+    typedef typename GET_PROP_TYPE(TypeTag, Grid) Grid;
+    typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
+    typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
+
+    typedef typename GET_PROP_TYPE(TypeTag, SolutionVector) SolutionVector;
+
+
+    typedef typename GET_PROP_TYPE(TypeTag, FVElementGeometry) FVElementGeometry;
+    typedef typename GridView::template Codim<0>::Entity Element;
+
+    //typedef LinearMaterial<Scalar> EffMaterialLaw;
+public:
+  
+    typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
+    typedef typename MaterialLaw::Params MaterialLawParams;
+    
+    RichardsNISpatialParams(const GridView &gridView)
+        : ParentType(gridView)
+    {
+        permeability_ = 1e-10;
+        porosity_ = 0.4;
+
+        // heat conductivity of granite
+        lambdaSolid_ = 2.8;
+	
+        // residual saturations
+
+        // residual saturations
+        materialParams_.setSwr(0.05);
+        materialParams_.setSnr(0.0);
+
+        // parameters for the Van Genuchten law
+        // alpha and n
+
+        materialParams_.setVgAlpha(0.0037);
+        materialParams_.setVgn(4.7);
+    }
+
+    ~RichardsNISpatialParams()
+    {}
+
+
+    /*!
+     * \brief Update the spatial parameters with the flow solution
+     *        after a timestep.
+     *
+     * \param globalSolution the global solution vector
+     */
+    void update(const SolutionVector &globalSolution)
+    {
+    }
+
+    /*!
+     * \brief Define the intrinsic permeability \f$\mathrm{[m^2]}\f$.
+     *
+     * \param element The current finite element
+     * \param fvGeometry The current finite volume geometry of the element
+     * \param scvIdx The index of the sub-control volume
+     */
+    const Scalar intrinsicPermeability(const Element &element,
+            const FVElementGeometry &fvGeometry,
+            const int scvIdx) const
+    {
+        return permeability_;
+    }
+
+    /*!
+     * \brief Define the porosity \f$\mathrm{[-]}\f$.
+     *
+     * \param element The finite element
+     * \param fvGeometry The finite volume geometry
+     * \param scvIdx The local index of the sub-control volume where
+     */
+    double porosity(const Element &element,
+            const FVElementGeometry &fvGeometry,
+            const int scvIdx) const
+    {
+        return porosity_;
+    }
+    
+        /*!
+     * \brief return the parameter object for the Brooks-Corey material law which depends on the position
+     *
+     * \param element The current finite element
+     * \param fvGeometry The current finite volume geometry of the element
+     * \param scvIdx The index of the sub-control volume
+     */
+    const MaterialLawParams& materialLawParams(const Element &element,
+            const FVElementGeometry &fvGeometry,
+            const int scvIdx) const
+    {
+        return materialParams_;
+    }
+
+
+    bool useTwoPointGradient(const Element &element,
+                             const int vertexI,
+                             const int vertexJ) const
+    {
+        return false;
+    }
+
+    /*!
+     * \brief Returns the heat capacity \f$[J/m^3 K]\f$ of the rock matrix.
+     *
+     * This is only required for non-isothermal models.
+     *
+     * \param element The finite element
+     * \param fvGeometry The finite volume geometry
+     * \param scvIdx The local index of the sub-control volume where
+     *                    the heat capacity needs to be defined
+     */
+    Scalar heatCapacity(const Element &element,
+                        const FVElementGeometry &fvGeometry,
+                        const int scvIdx) const
+    {
+        return
+            790 // specific heat capacity of granite [J / (kg K)]
+            * 2700 // density of granite [kg/m^3]
+            * (1 - porosity(element, fvGeometry, scvIdx));
+    }
+
+    /*!
+     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the porous material.
+     *
+     * \param element The finite element
+     * \param fvGeometry The finite volume geometry
+     * \param scvIdx The local index of the sub-control volume where
+     *                    the heat capacity needs to be defined
+     */
+    Scalar thermalConductivitySolid(const Element &element,
+                                    const FVElementGeometry &fvGeometry,
+                                    const int scvIdx) const
+    {
+        return lambdaSolid_;
+    }
+
+
+private:
+  
+    MaterialLawParams materialParams_;
+    Scalar permeability_;
+    Scalar porosity_;
+    Scalar lambdaSolid_;
+};
+
+}
+
+#endif
diff --git a/test/implicit/richards/test_boxrichardsniconduction.cc b/test/implicit/richards/test_boxrichardsniconduction.cc
new file mode 100644
index 0000000000000000000000000000000000000000..503b8defc7916fb6994c0b51bc37f7844c3dd72f
--- /dev/null
+++ b/test/implicit/richards/test_boxrichardsniconduction.cc
@@ -0,0 +1,57 @@
+// -*- 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 test for the RichardsNI box model
+ */
+#include "config.h"
+#include "richardsniconductionproblem.hh"
+#include <dumux/common/start.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ *        reading in parameters.
+ *
+ * \param progName  The name of the program, that was tried to be started.
+ * \param errorMsg  The error message that was issued by the start function.
+ *                  Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+    if (errorMsg.size() > 0) {
+        std::string errorMessageOut = "\nUsage: ";
+                    errorMessageOut += progName;
+                    errorMessageOut += " [options]\n";
+                    errorMessageOut += errorMsg;
+                    errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
+                                        "\t-TimeManager.TEnd      End of the simulation [s] \n"
+                                        "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+                                        "\t-Grid.File             Name of the file containing the grid \n"
+                                        "\t                       definition in DGF format\n";
+        std::cout << errorMessageOut
+                  << "\n";
+    }
+}
+
+int main(int argc, char** argv)
+{
+    typedef TTAG(RichardsNIConductionBoxProblem) ProblemTypeTag;
+    return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+}
diff --git a/test/implicit/richards/test_boxrichardsniconduction.input b/test/implicit/richards/test_boxrichardsniconduction.input
new file mode 100644
index 0000000000000000000000000000000000000000..72b6557d8aacce4b6972416cd8e7409ffb72eadb
--- /dev/null
+++ b/test/implicit/richards/test_boxrichardsniconduction.input
@@ -0,0 +1,12 @@
+[TimeManager]
+DtInitial = 1 # [s]
+TEnd = 1e5 # [s]
+MaxTimeStepSize = 1e10 # [s]
+
+[Grid]
+File = ./grids/test_richardsniconduction.dgf
+
+[Problem]
+Name = test_boxrichardsniconduction # name passed to the output routines
+OutputInterval = 5 # every 5th timestep an output file is written
+
diff --git a/test/implicit/richards/test_boxrichardsniconvection.cc b/test/implicit/richards/test_boxrichardsniconvection.cc
new file mode 100644
index 0000000000000000000000000000000000000000..5381c6ce0e5ee3fcdcdd5295b2520632561395c4
--- /dev/null
+++ b/test/implicit/richards/test_boxrichardsniconvection.cc
@@ -0,0 +1,57 @@
+// -*- 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 test for the RichardsNI box model
+ */
+#include "config.h"
+#include "richardsniconvectionproblem.hh"
+#include <dumux/common/start.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ *        reading in parameters.
+ *
+ * \param progName  The name of the program, that was tried to be started.
+ * \param errorMsg  The error message that was issued by the start function.
+ *                  Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+    if (errorMsg.size() > 0) {
+        std::string errorMessageOut = "\nUsage: ";
+                    errorMessageOut += progName;
+                    errorMessageOut += " [options]\n";
+                    errorMessageOut += errorMsg;
+                    errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
+                                        "\t-TimeManager.TEnd      End of the simulation [s] \n"
+                                        "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+                                        "\t-Grid.File             Name of the file containing the grid \n"
+                                        "\t                       definition in DGF format\n";
+        std::cout << errorMessageOut
+                  << "\n";
+    }
+}
+
+int main(int argc, char** argv)
+{
+    typedef TTAG(RichardsNIConvectionBoxProblem) ProblemTypeTag;
+    return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+}
diff --git a/test/implicit/richards/test_boxrichardsniconvection.input b/test/implicit/richards/test_boxrichardsniconvection.input
new file mode 100644
index 0000000000000000000000000000000000000000..47aa414db3df6d59a5dd2ad5e66a27c7eac386e1
--- /dev/null
+++ b/test/implicit/richards/test_boxrichardsniconvection.input
@@ -0,0 +1,12 @@
+[TimeManager]
+DtInitial = 1 # [s]
+TEnd = 3e4 # [s]
+MaxTimeStepSize = 1e3 # [s]
+
+[Grid]
+File = ./grids/test_richardsniconvection.dgf
+
+[Problem]
+Name = test_boxrichardsniconvection # name passed to the output routines
+OutputInterval = 5 # every 5th timestep an output file is written
+DarcyVelocity = 1e-4 # [m/s] inflow at the left boundary
diff --git a/test/implicit/richards/test_ccrichardsniconduction.cc b/test/implicit/richards/test_ccrichardsniconduction.cc
new file mode 100644
index 0000000000000000000000000000000000000000..18532cf19448bac3d107102e41c19b1e4adf9ff0
--- /dev/null
+++ b/test/implicit/richards/test_ccrichardsniconduction.cc
@@ -0,0 +1,57 @@
+// -*- 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 test for the RichardsNI cc model
+ */
+#include "config.h"
+#include "richardsniconductionproblem.hh"
+#include <dumux/common/start.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ *        reading in parameters.
+ *
+ * \param progName  The name of the program, that was tried to be started.
+ * \param errorMsg  The error message that was issued by the start function.
+ *                  Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+    if (errorMsg.size() > 0) {
+        std::string errorMessageOut = "\nUsage: ";
+                    errorMessageOut += progName;
+                    errorMessageOut += " [options]\n";
+                    errorMessageOut += errorMsg;
+                    errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
+                                        "\t-TimeManager.TEnd      End of the simulation [s] \n"
+                                        "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+                                        "\t-Grid.File             Name of the file containing the grid \n"
+                                        "\t                       definition in DGF format\n";
+        std::cout << errorMessageOut
+                  << "\n";
+    }
+}
+
+int main(int argc, char** argv)
+{
+    typedef TTAG(RichardsNIConductionCCProblem) ProblemTypeTag;
+    return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+}
diff --git a/test/implicit/richards/test_ccrichardsniconduction.input b/test/implicit/richards/test_ccrichardsniconduction.input
new file mode 100644
index 0000000000000000000000000000000000000000..096dc3a55959dd9108fb59d35e343d1044407dc2
--- /dev/null
+++ b/test/implicit/richards/test_ccrichardsniconduction.input
@@ -0,0 +1,11 @@
+[TimeManager]
+DtInitial = 1 # [s]
+TEnd = 1e5 # [s]
+MaxTimeStepSize = 1e10 # [s]
+
+[Grid]
+File = ./grids/test_richardsniconduction.dgf
+
+[Problem]
+Name = test_ccrichardsniconduction # name passed to the output routines
+OutputInterval = 5 # every 5th timestep an output file is written
diff --git a/test/implicit/richards/test_ccrichardsniconvection.cc b/test/implicit/richards/test_ccrichardsniconvection.cc
new file mode 100644
index 0000000000000000000000000000000000000000..bbd6edf30ae9dffa46876845fc033536da48650e
--- /dev/null
+++ b/test/implicit/richards/test_ccrichardsniconvection.cc
@@ -0,0 +1,57 @@
+// -*- 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 test for the RichardsNI cc model
+ */
+#include "config.h"
+#include "richardsniconvectionproblem.hh"
+#include <dumux/common/start.hh>
+
+/*!
+ * \brief Provides an interface for customizing error messages associated with
+ *        reading in parameters.
+ *
+ * \param progName  The name of the program, that was tried to be started.
+ * \param errorMsg  The error message that was issued by the start function.
+ *                  Comprises the thing that went wrong and a general help message.
+ */
+void usage(const char *progName, const std::string &errorMsg)
+{
+    if (errorMsg.size() > 0) {
+        std::string errorMessageOut = "\nUsage: ";
+                    errorMessageOut += progName;
+                    errorMessageOut += " [options]\n";
+                    errorMessageOut += errorMsg;
+                    errorMessageOut += "\n\nThe list of mandatory options for this program is:\n"
+                                        "\t-TimeManager.TEnd      End of the simulation [s] \n"
+                                        "\t-TimeManager.DtInitial Initial timestep size [s] \n"
+                                        "\t-Grid.File             Name of the file containing the grid \n"
+                                        "\t                       definition in DGF format\n";
+        std::cout << errorMessageOut
+                  << "\n";
+    }
+}
+
+int main(int argc, char** argv)
+{
+    typedef TTAG(RichardsNIConvectionCCProblem) ProblemTypeTag;
+    return Dumux::start<ProblemTypeTag>(argc, argv, usage);
+}
diff --git a/test/implicit/richards/test_ccrichardsniconvection.input b/test/implicit/richards/test_ccrichardsniconvection.input
new file mode 100644
index 0000000000000000000000000000000000000000..0bed72d6000594f8537b72c052285b487460a251
--- /dev/null
+++ b/test/implicit/richards/test_ccrichardsniconvection.input
@@ -0,0 +1,12 @@
+[TimeManager]
+DtInitial = 1 # [s]
+TEnd = 3e4 # [s]
+MaxTimeStepSize = 1e3 # [s]
+
+[Grid]
+File = ./grids/test_richardsniconvection.dgf
+
+[Problem]
+Name = test_ccrichardsniconvection # name passed to the output routines
+OutputInterval = 5 # every 5th timestep an output file is written
+DarcyVelocity = 1e-4 # [m/s] inflow at the left boundary
diff --git a/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cnispatialparams.hh b/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cnispatialparams.hh
index 0d9367fe19668a98efd0fffd17952164731ca626..580ae1f0a210f5700a304de86c04004051b9737c 100644
--- a/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cnispatialparams.hh
+++ b/test/multidomain/2cnistokes2p2cni/2cnistokes2p2cnispatialparams.hh
@@ -273,7 +273,7 @@ public:
     }
 
     /*!
-     * \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the solid
+     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the solid
      *
      * This is only required for non-isothermal models.
      *
diff --git a/test/multidomain/2cstokes2p2c/2cstokes2p2cspatialparams.hh b/test/multidomain/2cstokes2p2c/2cstokes2p2cspatialparams.hh
index 0bc423eb02fba01eb725173d50abdb0846a5da77..6d8e8105cbf1244cfb689571507cb996908b1f4f 100644
--- a/test/multidomain/2cstokes2p2c/2cstokes2p2cspatialparams.hh
+++ b/test/multidomain/2cstokes2p2c/2cstokes2p2cspatialparams.hh
@@ -279,7 +279,7 @@ public:
     }
 
     /*!
-     * \brief Returns the thermal conductivity \f$[W/m^2]\f$ of the solid
+     * \brief Returns the thermal conductivity \f$\mathrm{[W/(m K)]}\f$ of the solid
      *
      * This is only required for non-isothermal models.
      *
diff --git a/test/references/richardsniconductionbox-reference.vtu b/test/references/richardsniconductionbox-reference.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..f3fc0f5ab1a4f1495ff2ae85d5046fdce7875c7e
--- /dev/null
+++ b/test/references/richardsniconductionbox-reference.vtu
@@ -0,0 +1,779 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+  <UnstructuredGrid>
+    <Piece NumberOfCells="200" NumberOfPoints="402">
+      <PointData Scalars="Sn" Vectors="velocity">
+        <DataArray type="Float32" Name="Sn" 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 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="Sw" 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 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1 1 1 1 1 1 1
+          1 1 1 1 1 1
+        </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
+        </DataArray>
+        <DataArray type="Float32" Name="pw" 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
+        </DataArray>
+        <DataArray type="Float32" Name="pc" 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 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
+          996.544 996.726 996.544 996.726 996.906 996.906 997.084 997.084 997.239 997.239 997.385 997.385
+          997.526 997.526 997.66 997.66 997.787 997.787 997.906 997.906 998.007 998.007 998.094 998.094
+          998.175 998.175 998.248 998.248 998.315 998.315 998.376 998.376 998.431 998.431 998.48 998.48
+          998.524 998.524 998.562 998.562 998.597 998.597 998.627 998.627 998.654 998.654 998.673 998.673
+          998.69 998.69 998.704 998.704 998.716 998.716 998.727 998.727 998.736 998.736 998.744 998.744
+          998.751 998.751 998.757 998.757 998.762 998.762 998.766 998.766 998.77 998.77 998.773 998.773
+          998.775 998.775 998.777 998.777 998.779 998.779 998.781 998.781 998.782 998.782 998.783 998.783
+          998.784 998.784 998.785 998.785 998.786 998.786 998.786 998.786 998.787 998.787 998.787 998.787
+          998.788 998.788 998.788 998.788 998.788 998.788 998.788 998.788 998.788 998.788 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789
+        </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
+        </DataArray>
+        <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
+          1169.89 1152.46 1169.89 1152.46 1135.68 1135.68 1119.64 1119.64 1102.68 1102.68 1086.42 1086.42
+          1071.26 1071.26 1057.19 1057.19 1044.22 1044.22 1032.32 1032.32 1020.47 1020.47 1009.2 1009.2
+          999.031 999.031 989.899 989.899 981.729 981.729 974.446 974.446 967.977 967.977 962.249 962.249
+          957.194 957.194 952.746 952.746 948.845 948.845 945.431 945.431 942.453 942.453 939.496 939.496
+          936.9 936.9 934.653 934.653 932.713 932.713 931.041 931.041 929.604 929.604 928.371 928.371
+          927.314 927.314 926.411 926.411 925.64 925.64 924.983 924.983 924.424 924.424 923.949 923.949
+          923.546 923.546 923.205 923.205 922.917 922.917 922.673 922.673 922.467 922.467 922.293 922.293
+          922.147 922.147 922.024 922.024 921.921 921.921 921.834 921.834 921.762 921.762 921.701 921.701
+          921.65 921.65 921.607 921.607 921.571 921.571 921.541 921.541 921.516 921.516 921.495 921.495
+          921.478 921.478 921.464 921.464 921.452 921.452 921.442 921.442 921.433 921.433 921.426 921.426
+          921.42 921.42 921.416 921.416 921.412 921.412 921.408 921.408 921.406 921.406 921.403 921.403
+          921.401 921.401 921.4 921.4 921.398 921.398 921.397 921.397 921.396 921.396 921.396 921.396
+          921.395 921.395 921.395 921.395 921.394 921.394 921.394 921.394 921.394 921.394 921.393 921.393
+          921.393 921.393 921.393 921.393 921.393 921.393 921.393 921.393 921.393 921.393 921.393 921.393
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392
+        </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
+        </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
+          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
+          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
+          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
+          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
+          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
+          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
+          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
+          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
+          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
+          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
+          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
+          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
+          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
+          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
+          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
+          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">
+          300 299.327 300 299.327 298.659 298.659 298.003 298.003 297.362 297.362 296.742 296.742
+          296.147 296.147 295.58 295.58 295.043 295.043 294.539 294.539 294.068 294.068 293.632 293.632
+          293.23 293.23 292.863 292.863 292.528 292.528 292.225 292.225 291.952 291.952 291.707 291.707
+          291.488 291.488 291.294 291.294 291.122 291.122 290.97 290.97 290.837 290.837 290.721 290.721
+          290.619 290.619 290.531 290.531 290.454 290.454 290.388 290.388 290.33 290.33 290.281 290.281
+          290.239 290.239 290.203 290.203 290.172 290.172 290.145 290.145 290.123 290.123 290.104 290.104
+          290.087 290.087 290.073 290.073 290.062 290.062 290.052 290.052 290.044 290.044 290.037 290.037
+          290.031 290.031 290.026 290.026 290.021 290.021 290.018 290.018 290.015 290.015 290.013 290.013
+          290.01 290.01 290.009 290.009 290.007 290.007 290.006 290.006 290.005 290.005 290.004 290.004
+          290.003 290.003 290.003 290.003 290.002 290.002 290.002 290.002 290.002 290.002 290.001 290.001
+          290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290
+        </DataArray>
+        <DataArray type="Float32" Name="velocity" NumberOfComponents="3" format="ascii">
+          -1.04565e-09 -8.51209e-19 0 -1.03948e-09 6.28894e-19 0 -1.04565e-09 -8.51209e-19 0 -1.03948e-09 6.28894e-19 0
+          -1.02126e-09 4.13157e-19 0 -1.02126e-09 4.13157e-19 0 -9.91887e-10 4.68422e-18 0 -9.91887e-10 4.68422e-18 0
+          -9.54807e-10 4.61325e-18 0 -9.54807e-10 4.61325e-18 0 -9.11444e-10 2.56905e-18 0 -9.11444e-10 2.56905e-18 0
+          -8.61045e-10 4.09207e-18 0 -8.61045e-10 4.09207e-18 0 -8.04909e-10 3.65373e-18 0 -8.04909e-10 3.65373e-18 0
+          -7.44388e-10 1.3296e-18 0 -7.44388e-10 1.3296e-18 0 -6.8232e-10 -7.51112e-19 0 -6.8232e-10 -7.51112e-19 0
+          -6.23547e-10 -2.04186e-18 0 -6.23547e-10 -2.04186e-18 0 -5.67767e-10 -2.57001e-18 0 -5.67767e-10 -2.57001e-18 0
+          -5.12412e-10 -9.08613e-19 0 -5.12412e-10 -9.08613e-19 0 -4.58259e-10 1.62055e-18 0 -4.58259e-10 1.62055e-18 0
+          -4.05951e-10 1.78575e-18 0 -4.05951e-10 1.78575e-18 0 -3.55999e-10 0 0 -3.55999e-10 0 0
+          -3.08787e-10 1.76074e-18 0 -3.08787e-10 1.76074e-18 0 -2.64581e-10 1.75032e-19 0 -2.64581e-10 1.75032e-19 0
+          -2.23542e-10 3.48225e-18 0 -2.23542e-10 3.48225e-18 0 -1.85741e-10 3.29277e-18 0 -1.85741e-10 3.29277e-18 0
+          -1.51172e-10 2.41631e-18 0 -1.51172e-10 2.41631e-18 0 -1.19769e-10 2.40762e-18 0 -1.19769e-10 2.40762e-18 0
+          -9.38154e-11 3.42862e-18 0 -9.38154e-11 3.42862e-18 0 -7.31082e-11 1.53804e-18 0 -7.31082e-11 1.53804e-18 0
+          -5.48136e-11 3.40842e-18 0 -5.48136e-11 3.40842e-18 0 -3.85543e-11 1.53011e-18 0 -3.85543e-11 1.53011e-18 0
+          -2.4169e-11 3.90217e-18 0 -2.4169e-11 3.90217e-18 0 -1.14949e-11 1.69355e-19 0 -1.1495e-11 1.69355e-19 0
+          -3.72039e-13 -1.35275e-18 0 -3.72061e-13 -1.35275e-18 0 9.35382e-12 -1.6887e-19 0 9.35381e-12 -1.6887e-19 0
+          1.78291e-11 -1.01206e-18 0 1.78291e-11 -1.01206e-18 0 2.51914e-11 2.02216e-18 0 2.51914e-11 2.02216e-18 0
+          3.15677e-11 2.69397e-18 0 3.15677e-11 2.69397e-18 0 3.70748e-11 1.68253e-19 0 3.70748e-11 1.68253e-19 0
+          4.18188e-11 -1.51337e-18 0 4.18188e-11 -1.51337e-18 0 4.58955e-11 1.68065e-18 0 4.58955e-11 1.68065e-18 0
+          4.93907e-11 -1.67992e-19 0 4.93906e-11 -1.67992e-19 0 5.23809e-11 3.3586e-19 0 5.23809e-11 3.3586e-19 0
+          5.4934e-11 0 0 5.4934e-11 0 0 5.71097e-11 0 0 5.71097e-11 0 0
+          5.89605e-11 -2.85253e-18 0 5.89605e-11 -2.85253e-18 0 6.05323e-11 -3.18752e-18 0 6.05323e-11 -3.18752e-18 0
+          6.18651e-11 -3.01928e-18 0 6.18651e-11 -3.01928e-18 0 6.29934e-11 -3.18659e-18 0 6.29935e-11 -3.18659e-18 0
+          6.39474e-11 -1.00618e-18 0 6.39475e-11 -1.00618e-18 0 6.4753e-11 1.67681e-19 0 6.4753e-11 1.67681e-19 0
+          6.54324e-11 -1.84434e-18 0 6.54323e-11 -1.84434e-18 0 6.60046e-11 0 0 6.60046e-11 0 0
+          6.64861e-11 -1.50882e-18 0 6.6486e-11 -1.50882e-18 0 6.68908e-11 -1.34111e-18 0 6.68908e-11 -1.34111e-18 0
+          6.72306e-11 -2.68212e-18 0 6.72306e-11 -2.68212e-18 0 6.75157e-11 -2.01153e-18 0 6.75157e-11 -2.01153e-18 0
+          6.77547e-11 -1.34098e-18 0 6.77547e-11 -1.34098e-18 0 6.79548e-11 -2.01143e-18 0 6.79548e-11 -2.01143e-18 0
+          6.81223e-11 -3.01709e-18 0 6.81223e-11 -3.01709e-18 0 6.82624e-11 -5.69885e-18 0 6.82624e-11 -5.69885e-18 0
+          6.83795e-11 -5.69878e-18 0 6.83794e-11 -5.69878e-18 0 6.84772e-11 -3.18458e-18 0 6.84772e-11 -3.18458e-18 0
+          6.85587e-11 -5.36345e-18 0 6.85587e-11 -5.36345e-18 0 6.86268e-11 -4.86059e-18 0 6.86268e-11 -4.86059e-18 0
+          6.86835e-11 -4.69295e-18 0 6.86836e-11 -4.69295e-18 0 6.87308e-11 -1.50844e-18 0 6.87308e-11 -1.50844e-18 0
+          6.87702e-11 -4.02249e-18 0 6.87702e-11 -4.02249e-18 0 6.88029e-11 -1.17322e-18 0 6.88029e-11 -1.17322e-18 0
+          6.88302e-11 2.84925e-18 0 6.88302e-11 2.84925e-18 0 6.88529e-11 -5.02807e-19 0 6.88528e-11 -5.02807e-19 0
+          6.88717e-11 0 0 6.88717e-11 0 0 6.88873e-11 1.50841e-18 0 6.88874e-11 1.50841e-18 0
+          6.89004e-11 1.34081e-18 0 6.89004e-11 1.34081e-18 0 6.89112e-11 1.17321e-18 0 6.89112e-11 1.17321e-18 0
+          6.89202e-11 -1.50841e-18 0 6.89201e-11 -1.50841e-18 0 6.89276e-11 -1.67601e-18 0 6.89276e-11 -1.67601e-18 0
+          6.89338e-11 -1.34081e-18 0 6.89338e-11 -1.34081e-18 0 6.89389e-11 -1.67601e-19 0 6.89389e-11 -1.67601e-19 0
+          6.89431e-11 -1.1732e-18 0 6.89432e-11 -1.1732e-18 0 6.89467e-11 1.50841e-18 0 6.89467e-11 1.50841e-18 0
+          6.89496e-11 2.17881e-18 0 6.89496e-11 2.17881e-18 0 6.8952e-11 3.35201e-18 0 6.8952e-11 3.35201e-18 0
+          6.8954e-11 4.86041e-18 0 6.8954e-11 4.86041e-18 0 6.89557e-11 5.02801e-18 0 6.89557e-11 5.02801e-18 0
+          6.89571e-11 2.51401e-18 0 6.8957e-11 2.51401e-18 0 6.89582e-11 1.676e-18 0 6.89581e-11 1.676e-18 0
+          6.89591e-11 5.02801e-19 0 6.89591e-11 5.02801e-19 0 6.89599e-11 1.3408e-18 0 6.89599e-11 1.3408e-18 0
+          6.89606e-11 3.35201e-19 0 6.89606e-11 3.35201e-19 0 6.89611e-11 6.70401e-19 0 6.89611e-11 6.70401e-19 0
+          6.89615e-11 1.5084e-18 0 6.89615e-11 1.5084e-18 0 6.89619e-11 2.0112e-18 0 6.89619e-11 2.0112e-18 0
+          6.89622e-11 1.5084e-18 0 6.89622e-11 1.5084e-18 0 6.89625e-11 3.35201e-19 0 6.89625e-11 3.35201e-19 0
+          6.89626e-11 3.35201e-19 0 6.89626e-11 3.35201e-19 0 6.89628e-11 3.35201e-19 0 6.89628e-11 3.35201e-19 0
+          6.89629e-11 6.70401e-19 0 6.8963e-11 6.70401e-19 0 6.89631e-11 1.8436e-18 0 6.89631e-11 1.8436e-18 0
+          6.89632e-11 8.38001e-19 0 6.89632e-11 8.38001e-19 0 6.89633e-11 3.35201e-18 0 6.89633e-11 3.35201e-18 0
+          6.89633e-11 3.51961e-18 0 6.89633e-11 3.51961e-18 0 6.89634e-11 5.69841e-18 0 6.89634e-11 5.69841e-18 0
+          6.89634e-11 6.53641e-18 0 6.89634e-11 6.53641e-18 0 6.89635e-11 3.51961e-18 0 6.89635e-11 3.51961e-18 0
+          6.89635e-11 3.85481e-18 0 6.89635e-11 3.85481e-18 0 6.89635e-11 1.676e-19 0 6.89635e-11 1.676e-19 0
+          6.89636e-11 1.676e-19 0 6.89635e-11 1.676e-19 0 6.89636e-11 1.676e-19 0 6.89635e-11 1.676e-19 0
+          6.89636e-11 1.676e-19 0 6.89636e-11 1.676e-19 0 6.89636e-11 1.676e-19 0 6.89636e-11 1.676e-19 0
+          6.89636e-11 -5.02801e-19 0 6.89636e-11 -5.02801e-19 0 6.89636e-11 -1.676e-18 0 6.89636e-11 -1.676e-18 0
+          6.89636e-11 -3.35201e-19 0 6.89636e-11 -3.35201e-19 0 6.89636e-11 -3.35201e-19 0 6.89636e-11 -3.35201e-19 0
+          6.89636e-11 -2.6816e-18 0 6.89636e-11 -2.6816e-18 0 6.89636e-11 -1.676e-18 0 6.89637e-11 -1.676e-18 0
+          6.89636e-11 -1.676e-19 0 6.89636e-11 -1.676e-19 0 6.89636e-11 5.02801e-19 0 6.89636e-11 5.02801e-19 0
+          6.89636e-11 -1.3408e-18 0 6.89636e-11 -1.3408e-18 0 6.89636e-11 -1.676e-18 0 6.89636e-11 -1.676e-18 0
+          6.89637e-11 0 0 6.89637e-11 0 0 6.89636e-11 1.5084e-18 0 6.89636e-11 1.5084e-18 0
+          6.89636e-11 1.5084e-18 0 6.89636e-11 1.5084e-18 0 6.89636e-11 -3.35201e-19 0 6.89636e-11 -3.35201e-19 0
+          6.89636e-11 -1.5084e-18 0 6.89636e-11 -1.5084e-18 0 6.89636e-11 -1.3408e-18 0 6.89636e-11 -1.3408e-18 0
+          6.89636e-11 -1.676e-18 0 6.89636e-11 -1.676e-18 0 6.89636e-11 -1.676e-19 0 6.89636e-11 -1.676e-19 0
+          6.89636e-11 1.676e-19 0 6.89636e-11 1.676e-19 0 6.89636e-11 1.676e-18 0 6.89636e-11 1.676e-18 0
+          6.89636e-11 -1.676e-19 0 6.89636e-11 -1.676e-19 0 6.89636e-11 6.70401e-19 0 6.89636e-11 6.70401e-19 0
+          6.89636e-11 1.676e-19 0 6.89636e-11 1.676e-19 0 6.89636e-11 -3.35201e-19 0 6.89636e-11 -3.35201e-19 0
+          6.89636e-11 -1.676e-19 0 6.89636e-11 -1.676e-19 0 6.89636e-11 -1.8436e-18 0 6.89636e-11 -1.8436e-18 0
+          6.89636e-11 -3.35201e-19 0 6.89636e-11 -3.35201e-19 0 6.89636e-11 -1.676e-19 0 6.89636e-11 -1.676e-19 0
+          6.89636e-11 -1.3408e-18 0 6.89636e-11 -1.3408e-18 0 6.89636e-11 -1.5084e-18 0 6.89637e-11 -1.5084e-18 0
+          6.89636e-11 0 0 6.89636e-11 0 0 6.89636e-11 2.1788e-18 0 6.89636e-11 2.1788e-18 0
+          6.89636e-11 1.8436e-18 0 6.89636e-11 1.8436e-18 0 6.89636e-11 1.3408e-18 0 6.89636e-11 1.3408e-18 0
+          6.89637e-11 2.1788e-18 0 6.89636e-11 2.1788e-18 0 6.89636e-11 1.8436e-18 0 6.89636e-11 1.8436e-18 0
+          6.89636e-11 1.3408e-18 0 6.89636e-11 1.3408e-18 0 6.89636e-11 2.1788e-18 0 6.89636e-11 2.1788e-18 0
+          6.89636e-11 3.85481e-18 0 6.89636e-11 3.85481e-18 0 6.89637e-11 3.1844e-18 0 6.89636e-11 3.1844e-18 0
+          6.89637e-11 2.0112e-18 0 6.89636e-11 2.0112e-18 0 6.89636e-11 2.3464e-18 0 6.89636e-11 2.3464e-18 0
+          6.89636e-11 3.35201e-18 0 6.89636e-11 3.35201e-18 0 6.89636e-11 3.35201e-18 0 6.89636e-11 3.35201e-18 0
+          6.89636e-11 2.3464e-18 0 6.89636e-11 2.3464e-18 0 6.89636e-11 1.8436e-18 0 6.89637e-11 1.8436e-18 0
+          6.89636e-11 1.5084e-18 0 6.89636e-11 1.5084e-18 0 6.89636e-11 5.02801e-18 0 6.89636e-11 5.02801e-18 0
+          6.89636e-11 5.19561e-18 0 6.89636e-11 5.19561e-18 0 6.89636e-11 3.35201e-18 0 6.89636e-11 3.35201e-18 0
+          6.89636e-11 2.3464e-18 0 6.89637e-11 2.3464e-18 0 6.89636e-11 3.68721e-18 0 6.89636e-11 3.68721e-18 0
+          6.89636e-11 2.6816e-18 0 6.89636e-11 2.6816e-18 0 6.89637e-11 2.8492e-18 0 6.89636e-11 2.8492e-18 0
+          6.89636e-11 2.1788e-18 0 6.89636e-11 2.1788e-18 0 6.89636e-11 -1.676e-19 0 6.89636e-11 -1.676e-19 0
+          6.89636e-11 -1.676e-19 0 6.89636e-11 -1.676e-19 0 6.89636e-11 -1.3408e-18 0 6.89636e-11 -1.3408e-18 0
+          6.89636e-11 -1.8436e-18 0 6.89636e-11 -1.8436e-18 0 6.89636e-11 -3.35201e-19 0 6.89636e-11 -3.35201e-19 0
+          6.89637e-11 -3.35201e-18 0 6.89636e-11 -3.35201e-18 0 6.89636e-11 -2.514e-18 0 6.89636e-11 -2.514e-18 0
+          6.89636e-11 -1.5084e-18 0 6.89636e-11 -1.5084e-18 0 6.89636e-11 0 0 6.89636e-11 0 0
+          6.89636e-11 1.676e-19 0 6.89636e-11 1.676e-19 0 6.89636e-11 -3.35201e-19 0 6.89636e-11 -3.35201e-19 0
+          6.89636e-11 -1.5084e-18 0 6.89636e-11 -1.5084e-18 0 6.89636e-11 -1.1732e-18 0 6.89636e-11 -1.1732e-18 0
+          6.89636e-11 -1.676e-19 0 6.89636e-11 -1.676e-19 0 6.89636e-11 0 0 6.89636e-11 0 0
+          6.89636e-11 2.1788e-18 0 6.89636e-11 2.1788e-18 0 6.89636e-11 1.676e-18 0 6.89636e-11 1.676e-18 0
+          6.89636e-11 -1.676e-19 0 6.89636e-11 -1.676e-19 0 6.89636e-11 6.70401e-19 0 6.89636e-11 6.70401e-19 0
+          6.89636e-11 1.676e-19 0 6.89636e-11 1.676e-19 0 6.89636e-11 -3.35201e-19 0 6.89636e-11 -3.35201e-19 0
+          6.89636e-11 -1.676e-19 0 6.89636e-11 -1.676e-19 0 6.89636e-11 -1.8436e-18 0 6.89636e-11 -1.8436e-18 0
+          6.89636e-11 -3.35201e-19 0 6.89636e-11 -3.35201e-19 0 6.89636e-11 -1.676e-19 0 6.89636e-11 -1.676e-19 0
+          6.89636e-11 6.70401e-19 0 6.89636e-11 6.70401e-19 0 6.89636e-11 1.676e-19 0 6.89636e-11 1.676e-19 0
+          6.89636e-11 -3.35201e-19 0 6.89636e-11 -3.35201e-19 0 6.89637e-11 0 0 6.89636e-11 0 0
+          6.89636e-11 -6.70401e-19 0 6.89636e-11 -6.70401e-19 0 6.89636e-11 -1.676e-18 0 6.89636e-11 -1.676e-18 0
+          6.89637e-11 -1.676e-19 0 6.89636e-11 -1.676e-19 0 6.89636e-11 -1.1732e-18 0 6.89636e-11 -1.1732e-18 0
+          6.89636e-11 -3.35201e-19 0 6.89636e-11 -3.35201e-19 0 6.89636e-11 -1.676e-19 0 6.89636e-11 -1.676e-19 0
+          6.89636e-11 0 0 6.89636e-11 0 0 6.89636e-11 3.35201e-19 0 6.89636e-11 3.35201e-19 0
+          6.89636e-11 -1.676e-19 0 6.89636e-11 -1.676e-19 0 6.89636e-11 6.70401e-19 0 6.89636e-11 6.70401e-19 0
+          6.89636e-11 0 0 6.89636e-11 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="temperatureExact" NumberOfComponents="1" format="ascii">
+          300 299.373 300 299.373 298.749 298.749 298.133 298.133 297.529 297.529 296.939 296.939
+          296.368 296.368 295.817 295.817 295.289 295.289 294.787 294.787 294.313 294.313 293.866 293.866
+          293.449 293.449 293.062 293.062 292.705 292.705 292.378 292.378 292.079 292.079 291.809 291.809
+          291.566 291.566 291.348 291.348 291.155 291.155 290.984 290.984 290.834 290.834 290.703 290.703
+          290.589 290.589 290.491 290.491 290.407 290.407 290.336 290.336 290.275 290.275 290.225 290.225
+          290.182 290.182 290.147 290.147 290.118 290.118 290.094 290.094 290.075 290.075 290.059 290.059
+          290.046 290.046 290.036 290.036 290.028 290.028 290.021 290.021 290.016 290.016 290.013 290.013
+          290.009 290.009 290.007 290.007 290.005 290.005 290.004 290.004 290.003 290.003 290.002 290.002
+          290.002 290.002 290.001 290.001 290.001 290.001 290.001 290.001 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290
+        </DataArray>
+      </PointData>
+      <CellData Scalars="process rank">
+        <DataArray type="Float32" Name="process rank" 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
+        </DataArray>
+      </CellData>
+      <Points>
+        <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+          0 0 0 0.025 0 0 0 1 0 0.025 1 0
+          0.05 0 0 0.05 1 0 0.075 0 0 0.075 1 0
+          0.1 0 0 0.1 1 0 0.125 0 0 0.125 1 0
+          0.15 0 0 0.15 1 0 0.175 0 0 0.175 1 0
+          0.2 0 0 0.2 1 0 0.225 0 0 0.225 1 0
+          0.25 0 0 0.25 1 0 0.275 0 0 0.275 1 0
+          0.3 0 0 0.3 1 0 0.325 0 0 0.325 1 0
+          0.35 0 0 0.35 1 0 0.375 0 0 0.375 1 0
+          0.4 0 0 0.4 1 0 0.425 0 0 0.425 1 0
+          0.45 0 0 0.45 1 0 0.475 0 0 0.475 1 0
+          0.5 0 0 0.5 1 0 0.525 0 0 0.525 1 0
+          0.55 0 0 0.55 1 0 0.575 0 0 0.575 1 0
+          0.6 0 0 0.6 1 0 0.625 0 0 0.625 1 0
+          0.65 0 0 0.65 1 0 0.675 0 0 0.675 1 0
+          0.7 0 0 0.7 1 0 0.725 0 0 0.725 1 0
+          0.75 0 0 0.75 1 0 0.775 0 0 0.775 1 0
+          0.8 0 0 0.8 1 0 0.825 0 0 0.825 1 0
+          0.85 0 0 0.85 1 0 0.875 0 0 0.875 1 0
+          0.9 0 0 0.9 1 0 0.925 0 0 0.925 1 0
+          0.95 0 0 0.95 1 0 0.975 0 0 0.975 1 0
+          1 0 0 1 1 0 1.025 0 0 1.025 1 0
+          1.05 0 0 1.05 1 0 1.075 0 0 1.075 1 0
+          1.1 0 0 1.1 1 0 1.125 0 0 1.125 1 0
+          1.15 0 0 1.15 1 0 1.175 0 0 1.175 1 0
+          1.2 0 0 1.2 1 0 1.225 0 0 1.225 1 0
+          1.25 0 0 1.25 1 0 1.275 0 0 1.275 1 0
+          1.3 0 0 1.3 1 0 1.325 0 0 1.325 1 0
+          1.35 0 0 1.35 1 0 1.375 0 0 1.375 1 0
+          1.4 0 0 1.4 1 0 1.425 0 0 1.425 1 0
+          1.45 0 0 1.45 1 0 1.475 0 0 1.475 1 0
+          1.5 0 0 1.5 1 0 1.525 0 0 1.525 1 0
+          1.55 0 0 1.55 1 0 1.575 0 0 1.575 1 0
+          1.6 0 0 1.6 1 0 1.625 0 0 1.625 1 0
+          1.65 0 0 1.65 1 0 1.675 0 0 1.675 1 0
+          1.7 0 0 1.7 1 0 1.725 0 0 1.725 1 0
+          1.75 0 0 1.75 1 0 1.775 0 0 1.775 1 0
+          1.8 0 0 1.8 1 0 1.825 0 0 1.825 1 0
+          1.85 0 0 1.85 1 0 1.875 0 0 1.875 1 0
+          1.9 0 0 1.9 1 0 1.925 0 0 1.925 1 0
+          1.95 0 0 1.95 1 0 1.975 0 0 1.975 1 0
+          2 0 0 2 1 0 2.025 0 0 2.025 1 0
+          2.05 0 0 2.05 1 0 2.075 0 0 2.075 1 0
+          2.1 0 0 2.1 1 0 2.125 0 0 2.125 1 0
+          2.15 0 0 2.15 1 0 2.175 0 0 2.175 1 0
+          2.2 0 0 2.2 1 0 2.225 0 0 2.225 1 0
+          2.25 0 0 2.25 1 0 2.275 0 0 2.275 1 0
+          2.3 0 0 2.3 1 0 2.325 0 0 2.325 1 0
+          2.35 0 0 2.35 1 0 2.375 0 0 2.375 1 0
+          2.4 0 0 2.4 1 0 2.425 0 0 2.425 1 0
+          2.45 0 0 2.45 1 0 2.475 0 0 2.475 1 0
+          2.5 0 0 2.5 1 0 2.525 0 0 2.525 1 0
+          2.55 0 0 2.55 1 0 2.575 0 0 2.575 1 0
+          2.6 0 0 2.6 1 0 2.625 0 0 2.625 1 0
+          2.65 0 0 2.65 1 0 2.675 0 0 2.675 1 0
+          2.7 0 0 2.7 1 0 2.725 0 0 2.725 1 0
+          2.75 0 0 2.75 1 0 2.775 0 0 2.775 1 0
+          2.8 0 0 2.8 1 0 2.825 0 0 2.825 1 0
+          2.85 0 0 2.85 1 0 2.875 0 0 2.875 1 0
+          2.9 0 0 2.9 1 0 2.925 0 0 2.925 1 0
+          2.95 0 0 2.95 1 0 2.975 0 0 2.975 1 0
+          3 0 0 3 1 0 3.025 0 0 3.025 1 0
+          3.05 0 0 3.05 1 0 3.075 0 0 3.075 1 0
+          3.1 0 0 3.1 1 0 3.125 0 0 3.125 1 0
+          3.15 0 0 3.15 1 0 3.175 0 0 3.175 1 0
+          3.2 0 0 3.2 1 0 3.225 0 0 3.225 1 0
+          3.25 0 0 3.25 1 0 3.275 0 0 3.275 1 0
+          3.3 0 0 3.3 1 0 3.325 0 0 3.325 1 0
+          3.35 0 0 3.35 1 0 3.375 0 0 3.375 1 0
+          3.4 0 0 3.4 1 0 3.425 0 0 3.425 1 0
+          3.45 0 0 3.45 1 0 3.475 0 0 3.475 1 0
+          3.5 0 0 3.5 1 0 3.525 0 0 3.525 1 0
+          3.55 0 0 3.55 1 0 3.575 0 0 3.575 1 0
+          3.6 0 0 3.6 1 0 3.625 0 0 3.625 1 0
+          3.65 0 0 3.65 1 0 3.675 0 0 3.675 1 0
+          3.7 0 0 3.7 1 0 3.725 0 0 3.725 1 0
+          3.75 0 0 3.75 1 0 3.775 0 0 3.775 1 0
+          3.8 0 0 3.8 1 0 3.825 0 0 3.825 1 0
+          3.85 0 0 3.85 1 0 3.875 0 0 3.875 1 0
+          3.9 0 0 3.9 1 0 3.925 0 0 3.925 1 0
+          3.95 0 0 3.95 1 0 3.975 0 0 3.975 1 0
+          4 0 0 4 1 0 4.025 0 0 4.025 1 0
+          4.05 0 0 4.05 1 0 4.075 0 0 4.075 1 0
+          4.1 0 0 4.1 1 0 4.125 0 0 4.125 1 0
+          4.15 0 0 4.15 1 0 4.175 0 0 4.175 1 0
+          4.2 0 0 4.2 1 0 4.225 0 0 4.225 1 0
+          4.25 0 0 4.25 1 0 4.275 0 0 4.275 1 0
+          4.3 0 0 4.3 1 0 4.325 0 0 4.325 1 0
+          4.35 0 0 4.35 1 0 4.375 0 0 4.375 1 0
+          4.4 0 0 4.4 1 0 4.425 0 0 4.425 1 0
+          4.45 0 0 4.45 1 0 4.475 0 0 4.475 1 0
+          4.5 0 0 4.5 1 0 4.525 0 0 4.525 1 0
+          4.55 0 0 4.55 1 0 4.575 0 0 4.575 1 0
+          4.6 0 0 4.6 1 0 4.625 0 0 4.625 1 0
+          4.65 0 0 4.65 1 0 4.675 0 0 4.675 1 0
+          4.7 0 0 4.7 1 0 4.725 0 0 4.725 1 0
+          4.75 0 0 4.75 1 0 4.775 0 0 4.775 1 0
+          4.8 0 0 4.8 1 0 4.825 0 0 4.825 1 0
+          4.85 0 0 4.85 1 0 4.875 0 0 4.875 1 0
+          4.9 0 0 4.9 1 0 4.925 0 0 4.925 1 0
+          4.95 0 0 4.95 1 0 4.975 0 0 4.975 1 0
+          5 0 0 5 1 0
+        </DataArray>
+      </Points>
+      <Cells>
+        <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+          0 1 3 2 1 4 5 3 4 6 7 5
+          6 8 9 7 8 10 11 9 10 12 13 11
+          12 14 15 13 14 16 17 15 16 18 19 17
+          18 20 21 19 20 22 23 21 22 24 25 23
+          24 26 27 25 26 28 29 27 28 30 31 29
+          30 32 33 31 32 34 35 33 34 36 37 35
+          36 38 39 37 38 40 41 39 40 42 43 41
+          42 44 45 43 44 46 47 45 46 48 49 47
+          48 50 51 49 50 52 53 51 52 54 55 53
+          54 56 57 55 56 58 59 57 58 60 61 59
+          60 62 63 61 62 64 65 63 64 66 67 65
+          66 68 69 67 68 70 71 69 70 72 73 71
+          72 74 75 73 74 76 77 75 76 78 79 77
+          78 80 81 79 80 82 83 81 82 84 85 83
+          84 86 87 85 86 88 89 87 88 90 91 89
+          90 92 93 91 92 94 95 93 94 96 97 95
+          96 98 99 97 98 100 101 99 100 102 103 101
+          102 104 105 103 104 106 107 105 106 108 109 107
+          108 110 111 109 110 112 113 111 112 114 115 113
+          114 116 117 115 116 118 119 117 118 120 121 119
+          120 122 123 121 122 124 125 123 124 126 127 125
+          126 128 129 127 128 130 131 129 130 132 133 131
+          132 134 135 133 134 136 137 135 136 138 139 137
+          138 140 141 139 140 142 143 141 142 144 145 143
+          144 146 147 145 146 148 149 147 148 150 151 149
+          150 152 153 151 152 154 155 153 154 156 157 155
+          156 158 159 157 158 160 161 159 160 162 163 161
+          162 164 165 163 164 166 167 165 166 168 169 167
+          168 170 171 169 170 172 173 171 172 174 175 173
+          174 176 177 175 176 178 179 177 178 180 181 179
+          180 182 183 181 182 184 185 183 184 186 187 185
+          186 188 189 187 188 190 191 189 190 192 193 191
+          192 194 195 193 194 196 197 195 196 198 199 197
+          198 200 201 199 200 202 203 201 202 204 205 203
+          204 206 207 205 206 208 209 207 208 210 211 209
+          210 212 213 211 212 214 215 213 214 216 217 215
+          216 218 219 217 218 220 221 219 220 222 223 221
+          222 224 225 223 224 226 227 225 226 228 229 227
+          228 230 231 229 230 232 233 231 232 234 235 233
+          234 236 237 235 236 238 239 237 238 240 241 239
+          240 242 243 241 242 244 245 243 244 246 247 245
+          246 248 249 247 248 250 251 249 250 252 253 251
+          252 254 255 253 254 256 257 255 256 258 259 257
+          258 260 261 259 260 262 263 261 262 264 265 263
+          264 266 267 265 266 268 269 267 268 270 271 269
+          270 272 273 271 272 274 275 273 274 276 277 275
+          276 278 279 277 278 280 281 279 280 282 283 281
+          282 284 285 283 284 286 287 285 286 288 289 287
+          288 290 291 289 290 292 293 291 292 294 295 293
+          294 296 297 295 296 298 299 297 298 300 301 299
+          300 302 303 301 302 304 305 303 304 306 307 305
+          306 308 309 307 308 310 311 309 310 312 313 311
+          312 314 315 313 314 316 317 315 316 318 319 317
+          318 320 321 319 320 322 323 321 322 324 325 323
+          324 326 327 325 326 328 329 327 328 330 331 329
+          330 332 333 331 332 334 335 333 334 336 337 335
+          336 338 339 337 338 340 341 339 340 342 343 341
+          342 344 345 343 344 346 347 345 346 348 349 347
+          348 350 351 349 350 352 353 351 352 354 355 353
+          354 356 357 355 356 358 359 357 358 360 361 359
+          360 362 363 361 362 364 365 363 364 366 367 365
+          366 368 369 367 368 370 371 369 370 372 373 371
+          372 374 375 373 374 376 377 375 376 378 379 377
+          378 380 381 379 380 382 383 381 382 384 385 383
+          384 386 387 385 386 388 389 387 388 390 391 389
+          390 392 393 391 392 394 395 393 394 396 397 395
+          396 398 399 397 398 400 401 399
+        </DataArray>
+        <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+          4 8 12 16 20 24 28 32 36 40 44 48
+          52 56 60 64 68 72 76 80 84 88 92 96
+          100 104 108 112 116 120 124 128 132 136 140 144
+          148 152 156 160 164 168 172 176 180 184 188 192
+          196 200 204 208 212 216 220 224 228 232 236 240
+          244 248 252 256 260 264 268 272 276 280 284 288
+          292 296 300 304 308 312 316 320 324 328 332 336
+          340 344 348 352 356 360 364 368 372 376 380 384
+          388 392 396 400 404 408 412 416 420 424 428 432
+          436 440 444 448 452 456 460 464 468 472 476 480
+          484 488 492 496 500 504 508 512 516 520 524 528
+          532 536 540 544 548 552 556 560 564 568 572 576
+          580 584 588 592 596 600 604 608 612 616 620 624
+          628 632 636 640 644 648 652 656 660 664 668 672
+          676 680 684 688 692 696 700 704 708 712 716 720
+          724 728 732 736 740 744 748 752 756 760 764 768
+          772 776 780 784 788 792 796 800
+        </DataArray>
+        <DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9
+        </DataArray>
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+</VTKFile>
diff --git a/test/references/richardsniconductioncc-reference.vtu b/test/references/richardsniconductioncc-reference.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..a096b629e6dfcb912293fef2de3fa8fdbf087662
--- /dev/null
+++ b/test/references/richardsniconductioncc-reference.vtu
@@ -0,0 +1,522 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+  <UnstructuredGrid>
+    <Piece NumberOfCells="200" NumberOfPoints="402">
+      <CellData Scalars="Sn" Vectors="velocity">
+        <DataArray type="Float32" Name="Sn" 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
+        </DataArray>
+        <DataArray type="Float32" Name="Sw" 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
+        </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
+        </DataArray>
+        <DataArray type="Float32" Name="pw" 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
+        </DataArray>
+        <DataArray type="Float32" Name="pc" 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
+        </DataArray>
+        <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
+          996.635 996.816 996.995 997.163 997.312 997.456 997.593 997.724 997.847 997.96 998.051 998.135
+          998.212 998.282 998.346 998.404 998.456 998.502 998.544 998.58 998.613 998.641 998.664 998.682
+          998.697 998.711 998.722 998.732 998.74 998.748 998.754 998.76 998.764 998.768 998.771 998.774
+          998.777 998.779 998.78 998.782 998.783 998.784 998.785 998.786 998.786 998.787 998.787 998.788
+          998.788 998.788 998.788 998.788 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789 998.789 998.789
+        </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
+        </DataArray>
+        <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
+          1161.11 1144 1127.6 1111.26 1094.48 1078.77 1064.16 1050.64 1038.21 1026.63 1014.77 1004.04
+          994.394 985.742 978.014 971.136 965.035 959.642 954.889 950.713 947.055 943.859 940.897 938.104
+          935.684 933.594 931.791 930.24 928.909 927.768 926.792 925.959 925.25 924.646 924.134 923.699
+          923.331 923.02 922.758 922.536 922.35 922.193 922.061 921.951 921.858 921.781 921.716 921.662
+          921.616 921.578 921.547 921.521 921.499 921.481 921.466 921.453 921.443 921.434 921.427 921.421
+          921.416 921.412 921.408 921.405 921.403 921.401 921.4 921.398 921.397 921.396 921.396 921.395
+          921.395 921.394 921.394 921.393 921.393 921.393 921.393 921.393 921.393 921.393 921.393 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+          921.392 921.392 921.392 921.392 921.392 921.392 921.392 921.392
+        </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
+        </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
+          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
+          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
+          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
+          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
+          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
+          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
+          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
+          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">
+          299.663 298.993 298.331 297.682 297.052 296.444 295.863 295.31 294.79 294.302 293.849 293.43
+          293.045 292.693 292.374 292.085 291.826 291.594 291.388 291.204 291.043 290.9 290.775 290.666
+          290.571 290.489 290.417 290.356 290.303 290.257 290.218 290.185 290.156 290.132 290.111 290.093
+          290.079 290.066 290.055 290.046 290.039 290.033 290.027 290.023 290.019 290.016 290.013 290.011
+          290.009 290.008 290.006 290.005 290.004 290.004 290.003 290.002 290.002 290.002 290.001 290.001
+          290.001 290.001 290.001 290.001 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290
+        </DataArray>
+        <DataArray type="Float32" Name="velocity" NumberOfComponents="3" format="ascii">
+          -1.03797e-09 0 0 -1.02573e-09 0 0 -1.00183e-09 0 0 -9.69465e-10 0 0
+          -9.30241e-10 0 0 -8.83398e-10 0 0 -8.30174e-10 0 0 -7.7189e-10 0 0
+          -7.09893e-10 0 0 -6.50497e-10 0 0 -5.94981e-10 0 0 -5.39432e-10 0 0
+          -4.84687e-10 0 0 -4.31453e-10 0 0 -3.80308e-10 0 0 -3.31701e-10 0 0
+          -2.85956e-10 0 0 -2.4329e-10 0 0 -2.0382e-10 0 0 -1.6758e-10 0 0
+          -1.34542e-10 0 0 -1.06341e-10 0 0 -8.36012e-11 0 0 -6.4169e-11 0 0
+          -4.68528e-11 0 0 -3.14953e-11 0 0 -1.79351e-11 0 0 -6.01105e-12 0 0
+          4.43408e-12 0 0 1.35505e-11 0 0 2.14803e-11 0 0 2.83564e-11 0 0
+          3.4301e-11 0 0 3.94261e-11 0 0 4.38332e-11 0 0 4.76134e-11 0 0
+          5.08484e-11 0 0 5.36109e-11 0 0 5.59649e-11 0 0 5.7967e-11 0 0
+          5.96667e-11 0 0 6.11072e-11 0 0 6.23261e-11 0 0 6.33557e-11 0 0
+          6.42243e-11 0 0 6.4956e-11 0 0 6.55716e-11 0 0 6.60889e-11 0 0
+          6.6523e-11 0 0 6.6887e-11 0 0 6.71918e-11 0 0 6.74469e-11 0 0
+          6.766e-11 0 0 6.78381e-11 0 0 6.79866e-11 0 0 6.81105e-11 0 0
+          6.82136e-11 0 0 6.82995e-11 0 0 6.8371e-11 0 0 6.84304e-11 0 0
+          6.84798e-11 0 0 6.85208e-11 0 0 6.85548e-11 0 0 6.8583e-11 0 0
+          6.86064e-11 0 0 6.86258e-11 0 0 6.86419e-11 0 0 6.86551e-11 0 0
+          6.86662e-11 0 0 6.86753e-11 0 0 6.86828e-11 0 0 6.86891e-11 0 0
+          6.86942e-11 0 0 6.86984e-11 0 0 6.87019e-11 0 0 6.87048e-11 0 0
+          6.87072e-11 0 0 6.87092e-11 0 0 6.87109e-11 0 0 6.87121e-11 0 0
+          6.87132e-11 0 0 6.87142e-11 0 0 6.8715e-11 0 0 6.87156e-11 0 0
+          6.87161e-11 0 0 6.87165e-11 0 0 6.87169e-11 0 0 6.87171e-11 0 0
+          6.87174e-11 0 0 6.87176e-11 0 0 6.87177e-11 0 0 6.87179e-11 0 0
+          6.8718e-11 0 0 6.8718e-11 0 0 6.87181e-11 0 0 6.87182e-11 0 0
+          6.87182e-11 0 0 6.87182e-11 0 0 6.87183e-11 0 0 6.87183e-11 0 0
+          6.87184e-11 0 0 6.87184e-11 0 0 6.87184e-11 0 0 6.87184e-11 0 0
+          6.87184e-11 0 0 6.87184e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0
+          6.87185e-11 0 0 6.87184e-11 0 0 6.87185e-11 0 0 6.87184e-11 0 0
+          6.87184e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0
+          6.87184e-11 0 0 6.87184e-11 0 0 6.87184e-11 0 0 6.87184e-11 0 0
+          6.87184e-11 0 0 6.87184e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0
+          6.87185e-11 0 0 6.87185e-11 0 0 6.87184e-11 0 0 6.87185e-11 0 0
+          6.87184e-11 0 0 6.87184e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0
+          6.87185e-11 0 0 6.87185e-11 0 0 6.87184e-11 0 0 6.87185e-11 0 0
+          6.87185e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0
+          6.87184e-11 0 0 6.87184e-11 0 0 6.87184e-11 0 0 6.87184e-11 0 0
+          6.87184e-11 0 0 6.87185e-11 0 0 6.87184e-11 0 0 6.87184e-11 0 0
+          6.87184e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0
+          6.87185e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0
+          6.87185e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0
+          6.87185e-11 0 0 6.87185e-11 0 0 6.87184e-11 0 0 6.87185e-11 0 0
+          6.87185e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0
+          6.87185e-11 0 0 6.87184e-11 0 0 6.87184e-11 0 0 6.87184e-11 0 0
+          6.87185e-11 0 0 6.87185e-11 0 0 6.87184e-11 0 0 6.87184e-11 0 0
+          6.87185e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0
+          6.87185e-11 0 0 6.87184e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0
+          6.87185e-11 0 0 6.87185e-11 0 0 6.87184e-11 0 0 6.87184e-11 0 0
+          6.87185e-11 0 0 6.87184e-11 0 0 6.87184e-11 0 0 6.87184e-11 0 0
+          6.87184e-11 0 0 6.87185e-11 0 0 6.87184e-11 0 0 6.87184e-11 0 0
+          6.87184e-11 0 0 6.87184e-11 0 0 6.87185e-11 0 0 6.87185e-11 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="process rank" 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
+        </DataArray>
+        <DataArray type="Float32" Name="temperatureExact" NumberOfComponents="1" format="ascii">
+          299.686 299.06 298.44 297.83 297.232 296.651 296.09 295.55 295.035 294.547 294.086 293.654
+          293.252 292.88 292.538 292.225 291.941 291.684 291.454 291.249 291.067 290.906 290.766 290.644
+          290.538 290.448 290.37 290.304 290.249 290.202 290.164 290.132 290.105 290.084 290.066 290.052
+          290.041 290.032 290.024 290.019 290.014 290.011 290.008 290.006 290.005 290.003 290.003 290.002
+          290.001 290.001 290.001 290.001 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290
+        </DataArray>
+      </CellData>
+      <Points>
+        <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+          0 0 0 0.025 0 0 0 1 0 0.025 1 0
+          0.05 0 0 0.05 1 0 0.075 0 0 0.075 1 0
+          0.1 0 0 0.1 1 0 0.125 0 0 0.125 1 0
+          0.15 0 0 0.15 1 0 0.175 0 0 0.175 1 0
+          0.2 0 0 0.2 1 0 0.225 0 0 0.225 1 0
+          0.25 0 0 0.25 1 0 0.275 0 0 0.275 1 0
+          0.3 0 0 0.3 1 0 0.325 0 0 0.325 1 0
+          0.35 0 0 0.35 1 0 0.375 0 0 0.375 1 0
+          0.4 0 0 0.4 1 0 0.425 0 0 0.425 1 0
+          0.45 0 0 0.45 1 0 0.475 0 0 0.475 1 0
+          0.5 0 0 0.5 1 0 0.525 0 0 0.525 1 0
+          0.55 0 0 0.55 1 0 0.575 0 0 0.575 1 0
+          0.6 0 0 0.6 1 0 0.625 0 0 0.625 1 0
+          0.65 0 0 0.65 1 0 0.675 0 0 0.675 1 0
+          0.7 0 0 0.7 1 0 0.725 0 0 0.725 1 0
+          0.75 0 0 0.75 1 0 0.775 0 0 0.775 1 0
+          0.8 0 0 0.8 1 0 0.825 0 0 0.825 1 0
+          0.85 0 0 0.85 1 0 0.875 0 0 0.875 1 0
+          0.9 0 0 0.9 1 0 0.925 0 0 0.925 1 0
+          0.95 0 0 0.95 1 0 0.975 0 0 0.975 1 0
+          1 0 0 1 1 0 1.025 0 0 1.025 1 0
+          1.05 0 0 1.05 1 0 1.075 0 0 1.075 1 0
+          1.1 0 0 1.1 1 0 1.125 0 0 1.125 1 0
+          1.15 0 0 1.15 1 0 1.175 0 0 1.175 1 0
+          1.2 0 0 1.2 1 0 1.225 0 0 1.225 1 0
+          1.25 0 0 1.25 1 0 1.275 0 0 1.275 1 0
+          1.3 0 0 1.3 1 0 1.325 0 0 1.325 1 0
+          1.35 0 0 1.35 1 0 1.375 0 0 1.375 1 0
+          1.4 0 0 1.4 1 0 1.425 0 0 1.425 1 0
+          1.45 0 0 1.45 1 0 1.475 0 0 1.475 1 0
+          1.5 0 0 1.5 1 0 1.525 0 0 1.525 1 0
+          1.55 0 0 1.55 1 0 1.575 0 0 1.575 1 0
+          1.6 0 0 1.6 1 0 1.625 0 0 1.625 1 0
+          1.65 0 0 1.65 1 0 1.675 0 0 1.675 1 0
+          1.7 0 0 1.7 1 0 1.725 0 0 1.725 1 0
+          1.75 0 0 1.75 1 0 1.775 0 0 1.775 1 0
+          1.8 0 0 1.8 1 0 1.825 0 0 1.825 1 0
+          1.85 0 0 1.85 1 0 1.875 0 0 1.875 1 0
+          1.9 0 0 1.9 1 0 1.925 0 0 1.925 1 0
+          1.95 0 0 1.95 1 0 1.975 0 0 1.975 1 0
+          2 0 0 2 1 0 2.025 0 0 2.025 1 0
+          2.05 0 0 2.05 1 0 2.075 0 0 2.075 1 0
+          2.1 0 0 2.1 1 0 2.125 0 0 2.125 1 0
+          2.15 0 0 2.15 1 0 2.175 0 0 2.175 1 0
+          2.2 0 0 2.2 1 0 2.225 0 0 2.225 1 0
+          2.25 0 0 2.25 1 0 2.275 0 0 2.275 1 0
+          2.3 0 0 2.3 1 0 2.325 0 0 2.325 1 0
+          2.35 0 0 2.35 1 0 2.375 0 0 2.375 1 0
+          2.4 0 0 2.4 1 0 2.425 0 0 2.425 1 0
+          2.45 0 0 2.45 1 0 2.475 0 0 2.475 1 0
+          2.5 0 0 2.5 1 0 2.525 0 0 2.525 1 0
+          2.55 0 0 2.55 1 0 2.575 0 0 2.575 1 0
+          2.6 0 0 2.6 1 0 2.625 0 0 2.625 1 0
+          2.65 0 0 2.65 1 0 2.675 0 0 2.675 1 0
+          2.7 0 0 2.7 1 0 2.725 0 0 2.725 1 0
+          2.75 0 0 2.75 1 0 2.775 0 0 2.775 1 0
+          2.8 0 0 2.8 1 0 2.825 0 0 2.825 1 0
+          2.85 0 0 2.85 1 0 2.875 0 0 2.875 1 0
+          2.9 0 0 2.9 1 0 2.925 0 0 2.925 1 0
+          2.95 0 0 2.95 1 0 2.975 0 0 2.975 1 0
+          3 0 0 3 1 0 3.025 0 0 3.025 1 0
+          3.05 0 0 3.05 1 0 3.075 0 0 3.075 1 0
+          3.1 0 0 3.1 1 0 3.125 0 0 3.125 1 0
+          3.15 0 0 3.15 1 0 3.175 0 0 3.175 1 0
+          3.2 0 0 3.2 1 0 3.225 0 0 3.225 1 0
+          3.25 0 0 3.25 1 0 3.275 0 0 3.275 1 0
+          3.3 0 0 3.3 1 0 3.325 0 0 3.325 1 0
+          3.35 0 0 3.35 1 0 3.375 0 0 3.375 1 0
+          3.4 0 0 3.4 1 0 3.425 0 0 3.425 1 0
+          3.45 0 0 3.45 1 0 3.475 0 0 3.475 1 0
+          3.5 0 0 3.5 1 0 3.525 0 0 3.525 1 0
+          3.55 0 0 3.55 1 0 3.575 0 0 3.575 1 0
+          3.6 0 0 3.6 1 0 3.625 0 0 3.625 1 0
+          3.65 0 0 3.65 1 0 3.675 0 0 3.675 1 0
+          3.7 0 0 3.7 1 0 3.725 0 0 3.725 1 0
+          3.75 0 0 3.75 1 0 3.775 0 0 3.775 1 0
+          3.8 0 0 3.8 1 0 3.825 0 0 3.825 1 0
+          3.85 0 0 3.85 1 0 3.875 0 0 3.875 1 0
+          3.9 0 0 3.9 1 0 3.925 0 0 3.925 1 0
+          3.95 0 0 3.95 1 0 3.975 0 0 3.975 1 0
+          4 0 0 4 1 0 4.025 0 0 4.025 1 0
+          4.05 0 0 4.05 1 0 4.075 0 0 4.075 1 0
+          4.1 0 0 4.1 1 0 4.125 0 0 4.125 1 0
+          4.15 0 0 4.15 1 0 4.175 0 0 4.175 1 0
+          4.2 0 0 4.2 1 0 4.225 0 0 4.225 1 0
+          4.25 0 0 4.25 1 0 4.275 0 0 4.275 1 0
+          4.3 0 0 4.3 1 0 4.325 0 0 4.325 1 0
+          4.35 0 0 4.35 1 0 4.375 0 0 4.375 1 0
+          4.4 0 0 4.4 1 0 4.425 0 0 4.425 1 0
+          4.45 0 0 4.45 1 0 4.475 0 0 4.475 1 0
+          4.5 0 0 4.5 1 0 4.525 0 0 4.525 1 0
+          4.55 0 0 4.55 1 0 4.575 0 0 4.575 1 0
+          4.6 0 0 4.6 1 0 4.625 0 0 4.625 1 0
+          4.65 0 0 4.65 1 0 4.675 0 0 4.675 1 0
+          4.7 0 0 4.7 1 0 4.725 0 0 4.725 1 0
+          4.75 0 0 4.75 1 0 4.775 0 0 4.775 1 0
+          4.8 0 0 4.8 1 0 4.825 0 0 4.825 1 0
+          4.85 0 0 4.85 1 0 4.875 0 0 4.875 1 0
+          4.9 0 0 4.9 1 0 4.925 0 0 4.925 1 0
+          4.95 0 0 4.95 1 0 4.975 0 0 4.975 1 0
+          5 0 0 5 1 0
+        </DataArray>
+      </Points>
+      <Cells>
+        <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+          0 1 3 2 1 4 5 3 4 6 7 5
+          6 8 9 7 8 10 11 9 10 12 13 11
+          12 14 15 13 14 16 17 15 16 18 19 17
+          18 20 21 19 20 22 23 21 22 24 25 23
+          24 26 27 25 26 28 29 27 28 30 31 29
+          30 32 33 31 32 34 35 33 34 36 37 35
+          36 38 39 37 38 40 41 39 40 42 43 41
+          42 44 45 43 44 46 47 45 46 48 49 47
+          48 50 51 49 50 52 53 51 52 54 55 53
+          54 56 57 55 56 58 59 57 58 60 61 59
+          60 62 63 61 62 64 65 63 64 66 67 65
+          66 68 69 67 68 70 71 69 70 72 73 71
+          72 74 75 73 74 76 77 75 76 78 79 77
+          78 80 81 79 80 82 83 81 82 84 85 83
+          84 86 87 85 86 88 89 87 88 90 91 89
+          90 92 93 91 92 94 95 93 94 96 97 95
+          96 98 99 97 98 100 101 99 100 102 103 101
+          102 104 105 103 104 106 107 105 106 108 109 107
+          108 110 111 109 110 112 113 111 112 114 115 113
+          114 116 117 115 116 118 119 117 118 120 121 119
+          120 122 123 121 122 124 125 123 124 126 127 125
+          126 128 129 127 128 130 131 129 130 132 133 131
+          132 134 135 133 134 136 137 135 136 138 139 137
+          138 140 141 139 140 142 143 141 142 144 145 143
+          144 146 147 145 146 148 149 147 148 150 151 149
+          150 152 153 151 152 154 155 153 154 156 157 155
+          156 158 159 157 158 160 161 159 160 162 163 161
+          162 164 165 163 164 166 167 165 166 168 169 167
+          168 170 171 169 170 172 173 171 172 174 175 173
+          174 176 177 175 176 178 179 177 178 180 181 179
+          180 182 183 181 182 184 185 183 184 186 187 185
+          186 188 189 187 188 190 191 189 190 192 193 191
+          192 194 195 193 194 196 197 195 196 198 199 197
+          198 200 201 199 200 202 203 201 202 204 205 203
+          204 206 207 205 206 208 209 207 208 210 211 209
+          210 212 213 211 212 214 215 213 214 216 217 215
+          216 218 219 217 218 220 221 219 220 222 223 221
+          222 224 225 223 224 226 227 225 226 228 229 227
+          228 230 231 229 230 232 233 231 232 234 235 233
+          234 236 237 235 236 238 239 237 238 240 241 239
+          240 242 243 241 242 244 245 243 244 246 247 245
+          246 248 249 247 248 250 251 249 250 252 253 251
+          252 254 255 253 254 256 257 255 256 258 259 257
+          258 260 261 259 260 262 263 261 262 264 265 263
+          264 266 267 265 266 268 269 267 268 270 271 269
+          270 272 273 271 272 274 275 273 274 276 277 275
+          276 278 279 277 278 280 281 279 280 282 283 281
+          282 284 285 283 284 286 287 285 286 288 289 287
+          288 290 291 289 290 292 293 291 292 294 295 293
+          294 296 297 295 296 298 299 297 298 300 301 299
+          300 302 303 301 302 304 305 303 304 306 307 305
+          306 308 309 307 308 310 311 309 310 312 313 311
+          312 314 315 313 314 316 317 315 316 318 319 317
+          318 320 321 319 320 322 323 321 322 324 325 323
+          324 326 327 325 326 328 329 327 328 330 331 329
+          330 332 333 331 332 334 335 333 334 336 337 335
+          336 338 339 337 338 340 341 339 340 342 343 341
+          342 344 345 343 344 346 347 345 346 348 349 347
+          348 350 351 349 350 352 353 351 352 354 355 353
+          354 356 357 355 356 358 359 357 358 360 361 359
+          360 362 363 361 362 364 365 363 364 366 367 365
+          366 368 369 367 368 370 371 369 370 372 373 371
+          372 374 375 373 374 376 377 375 376 378 379 377
+          378 380 381 379 380 382 383 381 382 384 385 383
+          384 386 387 385 386 388 389 387 388 390 391 389
+          390 392 393 391 392 394 395 393 394 396 397 395
+          396 398 399 397 398 400 401 399
+        </DataArray>
+        <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+          4 8 12 16 20 24 28 32 36 40 44 48
+          52 56 60 64 68 72 76 80 84 88 92 96
+          100 104 108 112 116 120 124 128 132 136 140 144
+          148 152 156 160 164 168 172 176 180 184 188 192
+          196 200 204 208 212 216 220 224 228 232 236 240
+          244 248 252 256 260 264 268 272 276 280 284 288
+          292 296 300 304 308 312 316 320 324 328 332 336
+          340 344 348 352 356 360 364 368 372 376 380 384
+          388 392 396 400 404 408 412 416 420 424 428 432
+          436 440 444 448 452 456 460 464 468 472 476 480
+          484 488 492 496 500 504 508 512 516 520 524 528
+          532 536 540 544 548 552 556 560 564 568 572 576
+          580 584 588 592 596 600 604 608 612 616 620 624
+          628 632 636 640 644 648 652 656 660 664 668 672
+          676 680 684 688 692 696 700 704 708 712 716 720
+          724 728 732 736 740 744 748 752 756 760 764 768
+          772 776 780 784 788 792 796 800
+        </DataArray>
+        <DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9
+        </DataArray>
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+</VTKFile>
diff --git a/test/references/richardsniconvectionbox-reference.vtu b/test/references/richardsniconvectionbox-reference.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..1b068d0c803a5132b2a8a63bbc42019fe2babd63
--- /dev/null
+++ b/test/references/richardsniconvectionbox-reference.vtu
@@ -0,0 +1,349 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+  <UnstructuredGrid>
+    <Piece NumberOfCells="80" NumberOfPoints="162">
+      <PointData Scalars="Sn" Vectors="velocity">
+        <DataArray type="Float32" Name="Sn" 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
+        </DataArray>
+        <DataArray type="Float32" Name="Sw" 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
+        </DataArray>
+        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
+          121579 121315 121579 121315 121051 121051 120787 120787 120522 120522 120258 120258
+          119994 119994 119730 119730 119465 119465 119201 119201 118936 118936 118672 118672
+          118407 118407 118141 118141 117876 117876 117610 117610 117343 117343 117075 117075
+          116808 116808 116539 116539 116270 116270 116001 116001 115731 115731 115461 115461
+          115191 115191 114920 114920 114649 114649 114378 114378 114107 114107 113836 113836
+          113565 113565 113294 113294 113022 113022 112751 112751 112480 112480 112209 112209
+          111937 111937 111666 111666 111395 111395 111123 111123 110852 110852 110581 110581
+          110309 110309 110038 110038 109767 109767 109496 109496 109224 109224 108953 108953
+          108682 108682 108410 108410 108139 108139 107868 107868 107596 107596 107325 107325
+          107054 107054 106783 106783 106511 106511 106240 106240 105969 105969 105697 105697
+          105426 105426 105155 105155 104883 104883 104612 104612 104341 104341 104070 104070
+          103798 103798 103527 103527 103256 103256 102984 102984 102713 102713 102442 102442
+          102170 102170 101899 101899 101628 101628 101357 101357 101085 101085 100814 100814
+          100543 100543 100271 100271 100000 100000
+        </DataArray>
+        <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
+          121579 121315 121579 121315 121051 121051 120787 120787 120522 120522 120258 120258
+          119994 119994 119730 119730 119465 119465 119201 119201 118936 118936 118672 118672
+          118407 118407 118141 118141 117876 117876 117610 117610 117343 117343 117075 117075
+          116808 116808 116539 116539 116270 116270 116001 116001 115731 115731 115461 115461
+          115191 115191 114920 114920 114649 114649 114378 114378 114107 114107 113836 113836
+          113565 113565 113294 113294 113022 113022 112751 112751 112480 112480 112209 112209
+          111937 111937 111666 111666 111395 111395 111123 111123 110852 110852 110581 110581
+          110309 110309 110038 110038 109767 109767 109496 109496 109224 109224 108953 108953
+          108682 108682 108410 108410 108139 108139 107868 107868 107596 107596 107325 107325
+          107054 107054 106783 106783 106511 106511 106240 106240 105969 105969 105697 105697
+          105426 105426 105155 105155 104883 104883 104612 104612 104341 104341 104070 104070
+          103798 103798 103527 103527 103256 103256 102984 102984 102713 102713 102442 102442
+          102170 102170 101899 101899 101628 101628 101357 101357 101085 101085 100814 100814
+          100543 100543 100271 100271 100000 100000
+        </DataArray>
+        <DataArray type="Float32" Name="pc" 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
+        </DataArray>
+        <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
+          998.631 998.631 998.631 998.631 998.631 998.631 998.631 998.631 998.631 998.631 998.631 998.631
+          998.632 998.632 998.633 998.633 998.635 998.635 998.639 998.639 998.645 998.645 998.653 998.653
+          998.663 998.663 998.673 998.673 998.684 998.684 998.696 998.696 998.708 998.708 998.721 998.721
+          998.733 998.733 998.744 998.744 998.754 998.754 998.763 998.763 998.77 998.77 998.776 998.776
+          998.781 998.781 998.785 998.785 998.788 998.788 998.79 998.79 998.792 998.792 998.793 998.793
+          998.793 998.793 998.794 998.794 998.794 998.794 998.794 998.794 998.794 998.794 998.794 998.794
+          998.794 998.794 998.794 998.794 998.794 998.794 998.794 998.794 998.794 998.794 998.794 998.794
+          998.794 998.794 998.794 998.794 998.794 998.794 998.794 998.794 998.793 998.793 998.793 998.793
+          998.793 998.793 998.793 998.793 998.793 998.793 998.793 998.793 998.793 998.793 998.792 998.792
+          998.792 998.792 998.792 998.792 998.792 998.792 998.792 998.792 998.792 998.792 998.792 998.792
+          998.792 998.792 998.792 998.792 998.791 998.791 998.791 998.791 998.791 998.791 998.791 998.791
+          998.791 998.791 998.791 998.791 998.791 998.791 998.79 998.79 998.79 998.79 998.79 998.79
+          998.79 998.79 998.79 998.79 998.79 998.79 998.79 998.79 998.79 998.79 998.789 998.789
+          998.789 998.789 998.789 998.789 998.789 998.789
+        </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
+        </DataArray>
+        <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
+          946.105 946.107 946.105 946.107 946.107 946.107 946.104 946.104 946.093 946.093 946.059 946.059
+          945.981 945.981 945.825 945.825 945.548 945.548 945.104 945.104 944.452 944.452 943.563 943.563
+          942.432 942.432 940.895 940.895 939.113 939.113 937.185 937.185 935.189 935.189 933.206 933.206
+          931.309 931.309 929.553 929.553 927.981 927.981 926.613 926.613 925.456 925.456 924.501 924.501
+          923.731 923.731 923.126 923.126 922.659 922.659 922.307 922.307 922.046 922.046 921.856 921.856
+          921.72 921.72 921.624 921.624 921.558 921.558 921.512 921.512 921.482 921.482 921.462 921.462
+          921.448 921.448 921.44 921.44 921.434 921.434 921.43 921.43 921.428 921.428 921.427 921.427
+          921.426 921.426 921.425 921.425 921.425 921.425 921.424 921.424 921.424 921.424 921.424 921.424
+          921.424 921.424 921.423 921.423 921.423 921.423 921.423 921.423 921.423 921.423 921.423 921.423
+          921.423 921.423 921.422 921.422 921.422 921.422 921.422 921.422 921.422 921.422 921.422 921.422
+          921.422 921.422 921.422 921.422 921.421 921.421 921.421 921.421 921.421 921.421 921.421 921.421
+          921.421 921.421 921.421 921.421 921.421 921.421 921.42 921.42 921.42 921.42 921.42 921.42
+          921.42 921.42 921.42 921.42 921.42 921.42 921.419 921.419 921.419 921.419 921.419 921.419
+          921.419 921.419 921.419 921.419 921.392 921.392
+        </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
+        </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
+          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
+          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
+          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
+          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
+          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
+          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">
+          291 291 291 291 291 291 291 291 291 291 290.998 290.998
+          290.995 290.995 290.988 290.988 290.975 290.975 290.956 290.956 290.926 290.926 290.887 290.887
+          290.836 290.836 290.775 290.775 290.706 290.706 290.63 290.63 290.552 290.552 290.473 290.473
+          290.398 290.398 290.328 290.328 290.265 290.265 290.211 290.211 290.164 290.164 290.126 290.126
+          290.094 290.094 290.07 290.07 290.051 290.051 290.037 290.037 290.026 290.026 290.019 290.019
+          290.013 290.013 290.009 290.009 290.007 290.007 290.005 290.005 290.003 290.003 290.003 290.003
+          290.002 290.002 290.002 290.002 290.002 290.002 290.001 290.001 290.001 290.001 290.001 290.001
+          290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001
+          290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001
+          290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001
+          290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001
+          290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001
+          290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001
+          290.001 290.001 290.001 290.001 290 290
+        </DataArray>
+        <DataArray type="Float32" Name="velocity" NumberOfComponents="3" format="ascii">
+          0.0001 -3.44191e-19 0 0.0001 -5.16287e-19 0 0.0001 -3.44191e-19 0 0.0001 -5.16287e-19 0
+          0.0001 0 0 0.0001 0 0 0.0001 -3.44191e-19 0 0.0001 -3.44191e-19 0
+          0.0001 3.44187e-19 0 0.0001 3.44187e-19 0 0.0001 0 0 0.0001 0 0
+          0.0001 -1.54866e-18 0 0.0001 -1.54866e-18 0 0.0001 0 0 0.0001 0 0
+          9.99999e-05 -1.71994e-19 0 9.99999e-05 -1.71994e-19 0 9.99998e-05 1.71913e-19 0 9.99998e-05 1.71913e-19 0
+          9.99996e-05 1.54615e-18 0 9.99996e-05 1.54615e-18 0 9.99993e-05 -1.71633e-19 0 9.99993e-05 -1.71633e-19 0
+          9.99989e-05 1.71427e-19 0 9.99989e-05 1.71427e-19 0 9.99984e-05 2.05377e-18 0 9.99984e-05 2.05377e-18 0
+          9.9998e-05 1.36659e-18 0 9.9998e-05 1.36659e-18 0 9.99975e-05 1.36378e-18 0 9.99975e-05 1.36378e-18 0
+          9.9997e-05 -3.4022e-19 0 9.9997e-05 -3.4022e-19 0 9.99964e-05 1.69749e-19 0 9.99964e-05 1.69749e-19 0
+          9.99959e-05 -1.69404e-19 0 9.99959e-05 -1.69404e-19 0 9.99954e-05 1.18359e-18 0 9.99954e-05 1.18359e-18 0
+          9.9995e-05 -1.01279e-18 0 9.9995e-05 -1.01279e-18 0 9.99946e-05 3.371e-19 0 9.99946e-05 3.371e-19 0
+          9.99942e-05 6.73357e-19 0 9.99942e-05 6.73357e-19 0 9.99939e-05 1.68166e-18 0 9.99939e-05 1.68166e-18 0
+          9.99937e-05 -1.51223e-18 0 9.99937e-05 -1.51223e-18 0 9.99935e-05 -2.01499e-18 0 9.99935e-05 -2.01499e-18 0
+          9.99934e-05 1.67831e-19 0 9.99934e-05 1.67831e-19 0 9.99933e-05 -3.35533e-19 0 9.99933e-05 -3.35533e-19 0
+          9.99932e-05 -3.35438e-19 0 9.99932e-05 -3.35438e-19 0 9.99932e-05 -5.03054e-19 0 9.99932e-05 -5.03054e-19 0
+          9.99931e-05 -3.3532e-19 0 9.99931e-05 -3.3532e-19 0 9.99931e-05 -6.7057e-19 0 9.99931e-05 -6.7057e-19 0
+          9.99931e-05 -8.38152e-19 0 9.99931e-05 -8.38152e-19 0 9.99931e-05 -1.67622e-19 0 9.99931e-05 -1.67622e-19 0
+          9.99931e-05 -1.67617e-19 0 9.99931e-05 -1.67617e-19 0 9.99931e-05 -3.35226e-19 0 9.99931e-05 -3.35226e-19 0
+          9.99931e-05 5.02831e-19 0 9.99931e-05 5.02831e-19 0 9.99931e-05 1.67609e-19 0 9.99931e-05 1.67609e-19 0
+          9.99931e-05 -3.35216e-19 0 9.99931e-05 -3.35216e-19 0 9.99932e-05 1.67607e-19 0 9.99932e-05 1.67607e-19 0
+          9.99932e-05 -8.38034e-19 0 9.99932e-05 -8.38034e-19 0 9.99932e-05 -5.0282e-19 0 9.99932e-05 -5.0282e-19 0
+          9.99932e-05 1.00564e-18 0 9.99932e-05 1.00564e-18 0 9.99932e-05 -1.00564e-18 0 9.99932e-05 -1.00564e-18 0
+          9.99932e-05 1.67606e-19 0 9.99932e-05 1.67606e-19 0 9.99932e-05 -8.3803e-19 0 9.99932e-05 -8.3803e-19 0
+          9.99932e-05 0 0 9.99932e-05 0 0 9.99933e-05 1.84367e-18 0 9.99933e-05 1.84367e-18 0
+          9.99933e-05 -1.17324e-18 0 9.99933e-05 -1.17324e-18 0 9.99933e-05 1.67606e-18 0 9.99933e-05 1.67606e-18 0
+          9.99933e-05 3.35212e-19 0 9.99933e-05 3.35212e-19 0 9.99933e-05 1.67606e-19 0 9.99933e-05 1.67606e-19 0
+          9.99933e-05 1.67606e-19 0 9.99933e-05 1.67606e-19 0 9.99933e-05 -1.67606e-19 0 9.99933e-05 -1.67606e-19 0
+          9.99933e-05 5.02817e-19 0 9.99933e-05 5.02817e-19 0 9.99934e-05 -5.02817e-19 0 9.99934e-05 -5.02817e-19 0
+          9.99934e-05 6.70423e-19 0 9.99934e-05 6.70423e-19 0 9.99934e-05 -6.70423e-19 0 9.99934e-05 -6.70423e-19 0
+          9.99934e-05 1.34085e-18 0 9.99934e-05 1.34085e-18 0 9.99934e-05 0 0 9.99934e-05 0 0
+          9.99934e-05 1.67606e-18 0 9.99934e-05 1.67606e-18 0 9.99934e-05 1.67606e-19 0 9.99934e-05 1.67606e-19 0
+          9.99934e-05 -3.35211e-19 0 9.99934e-05 -3.35211e-19 0 9.99935e-05 -1.50845e-18 0 9.99935e-05 -1.50845e-18 0
+          9.99935e-05 1.67606e-19 0 9.99935e-05 1.67606e-19 0 9.99935e-05 6.70422e-19 0 9.99935e-05 6.70422e-19 0
+          9.99935e-05 -1.67605e-19 0 9.99935e-05 -1.67605e-19 0 9.99935e-05 1.50845e-18 0 9.99935e-05 1.50845e-18 0
+          9.99935e-05 -1.17324e-18 0 9.99935e-05 -1.17324e-18 0 9.99935e-05 -1.67605e-19 0 9.99935e-05 -1.67605e-19 0
+          9.99936e-05 -1.67605e-18 0 9.99936e-05 -1.67605e-18 0 9.99936e-05 1.67605e-19 0 9.99936e-05 1.67605e-19 0
+          9.99936e-05 -5.02816e-19 0 9.99936e-05 -5.02816e-19 0 9.99936e-05 -1.67605e-18 0 9.99936e-05 -1.67605e-18 0
+          9.99936e-05 -6.70421e-19 0 9.99936e-05 -6.70421e-19 0 9.99936e-05 1.00563e-18 0 9.99936e-05 1.00563e-18 0
+          9.99936e-05 -5.02816e-19 0 9.99936e-05 -5.02816e-19 0 9.99936e-05 0 0 9.99936e-05 0 0
+          9.99937e-05 0 0 9.99937e-05 0 0 9.99937e-05 0 0 9.99937e-05 0 0
+          9.99937e-05 0 0 9.99937e-05 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="temperatureExact" NumberOfComponents="1" format="ascii">
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 291 291 291 291 291 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290
+        </DataArray>
+      </PointData>
+      <CellData Scalars="process rank">
+        <DataArray type="Float32" Name="process rank" 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
+        </DataArray>
+      </CellData>
+      <Points>
+        <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+          0 0 0 0.25 0 0 0 1 0 0.25 1 0
+          0.5 0 0 0.5 1 0 0.75 0 0 0.75 1 0
+          1 0 0 1 1 0 1.25 0 0 1.25 1 0
+          1.5 0 0 1.5 1 0 1.75 0 0 1.75 1 0
+          2 0 0 2 1 0 2.25 0 0 2.25 1 0
+          2.5 0 0 2.5 1 0 2.75 0 0 2.75 1 0
+          3 0 0 3 1 0 3.25 0 0 3.25 1 0
+          3.5 0 0 3.5 1 0 3.75 0 0 3.75 1 0
+          4 0 0 4 1 0 4.25 0 0 4.25 1 0
+          4.5 0 0 4.5 1 0 4.75 0 0 4.75 1 0
+          5 0 0 5 1 0 5.25 0 0 5.25 1 0
+          5.5 0 0 5.5 1 0 5.75 0 0 5.75 1 0
+          6 0 0 6 1 0 6.25 0 0 6.25 1 0
+          6.5 0 0 6.5 1 0 6.75 0 0 6.75 1 0
+          7 0 0 7 1 0 7.25 0 0 7.25 1 0
+          7.5 0 0 7.5 1 0 7.75 0 0 7.75 1 0
+          8 0 0 8 1 0 8.25 0 0 8.25 1 0
+          8.5 0 0 8.5 1 0 8.75 0 0 8.75 1 0
+          9 0 0 9 1 0 9.25 0 0 9.25 1 0
+          9.5 0 0 9.5 1 0 9.75 0 0 9.75 1 0
+          10 0 0 10 1 0 10.25 0 0 10.25 1 0
+          10.5 0 0 10.5 1 0 10.75 0 0 10.75 1 0
+          11 0 0 11 1 0 11.25 0 0 11.25 1 0
+          11.5 0 0 11.5 1 0 11.75 0 0 11.75 1 0
+          12 0 0 12 1 0 12.25 0 0 12.25 1 0
+          12.5 0 0 12.5 1 0 12.75 0 0 12.75 1 0
+          13 0 0 13 1 0 13.25 0 0 13.25 1 0
+          13.5 0 0 13.5 1 0 13.75 0 0 13.75 1 0
+          14 0 0 14 1 0 14.25 0 0 14.25 1 0
+          14.5 0 0 14.5 1 0 14.75 0 0 14.75 1 0
+          15 0 0 15 1 0 15.25 0 0 15.25 1 0
+          15.5 0 0 15.5 1 0 15.75 0 0 15.75 1 0
+          16 0 0 16 1 0 16.25 0 0 16.25 1 0
+          16.5 0 0 16.5 1 0 16.75 0 0 16.75 1 0
+          17 0 0 17 1 0 17.25 0 0 17.25 1 0
+          17.5 0 0 17.5 1 0 17.75 0 0 17.75 1 0
+          18 0 0 18 1 0 18.25 0 0 18.25 1 0
+          18.5 0 0 18.5 1 0 18.75 0 0 18.75 1 0
+          19 0 0 19 1 0 19.25 0 0 19.25 1 0
+          19.5 0 0 19.5 1 0 19.75 0 0 19.75 1 0
+          20 0 0 20 1 0
+        </DataArray>
+      </Points>
+      <Cells>
+        <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+          0 1 3 2 1 4 5 3 4 6 7 5
+          6 8 9 7 8 10 11 9 10 12 13 11
+          12 14 15 13 14 16 17 15 16 18 19 17
+          18 20 21 19 20 22 23 21 22 24 25 23
+          24 26 27 25 26 28 29 27 28 30 31 29
+          30 32 33 31 32 34 35 33 34 36 37 35
+          36 38 39 37 38 40 41 39 40 42 43 41
+          42 44 45 43 44 46 47 45 46 48 49 47
+          48 50 51 49 50 52 53 51 52 54 55 53
+          54 56 57 55 56 58 59 57 58 60 61 59
+          60 62 63 61 62 64 65 63 64 66 67 65
+          66 68 69 67 68 70 71 69 70 72 73 71
+          72 74 75 73 74 76 77 75 76 78 79 77
+          78 80 81 79 80 82 83 81 82 84 85 83
+          84 86 87 85 86 88 89 87 88 90 91 89
+          90 92 93 91 92 94 95 93 94 96 97 95
+          96 98 99 97 98 100 101 99 100 102 103 101
+          102 104 105 103 104 106 107 105 106 108 109 107
+          108 110 111 109 110 112 113 111 112 114 115 113
+          114 116 117 115 116 118 119 117 118 120 121 119
+          120 122 123 121 122 124 125 123 124 126 127 125
+          126 128 129 127 128 130 131 129 130 132 133 131
+          132 134 135 133 134 136 137 135 136 138 139 137
+          138 140 141 139 140 142 143 141 142 144 145 143
+          144 146 147 145 146 148 149 147 148 150 151 149
+          150 152 153 151 152 154 155 153 154 156 157 155
+          156 158 159 157 158 160 161 159
+        </DataArray>
+        <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+          4 8 12 16 20 24 28 32 36 40 44 48
+          52 56 60 64 68 72 76 80 84 88 92 96
+          100 104 108 112 116 120 124 128 132 136 140 144
+          148 152 156 160 164 168 172 176 180 184 188 192
+          196 200 204 208 212 216 220 224 228 232 236 240
+          244 248 252 256 260 264 268 272 276 280 284 288
+          292 296 300 304 308 312 316 320
+        </DataArray>
+        <DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9
+        </DataArray>
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+</VTKFile>
diff --git a/test/references/richardsniconvectioncc-reference.vtu b/test/references/richardsniconvectioncc-reference.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..79088c9310d27ceb3475dc3b4c67bed4d3ec5769
--- /dev/null
+++ b/test/references/richardsniconvectioncc-reference.vtu
@@ -0,0 +1,242 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
+  <UnstructuredGrid>
+    <Piece NumberOfCells="80" NumberOfPoints="162">
+      <CellData Scalars="Sn" Vectors="velocity">
+        <DataArray type="Float32" Name="Sn" 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
+        </DataArray>
+        <DataArray type="Float32" Name="Sw" 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
+        </DataArray>
+        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
+          121447 121183 120919 120654 120390 120126 119862 119597 119333 119069 118804 118539
+          118274 118009 117743 117476 117209 116941 116673 116405 116136 115866 115596 115326
+          115055 114784 114514 114243 113972 113700 113429 113158 112887 112615 112344 112073
+          111802 111530 111259 110988 110716 110445 110174 109903 109631 109360 109089 108817
+          108546 108275 108003 107732 107461 107190 106918 106647 106376 106104 105833 105562
+          105290 105019 104748 104476 104205 103934 103663 103391 103120 102849 102577 102306
+          102035 101763 101492 101221 100950 100678 100407 100136
+        </DataArray>
+        <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
+          121447 121183 120919 120654 120390 120126 119862 119597 119333 119069 118804 118539
+          118274 118009 117743 117476 117209 116941 116673 116405 116136 115866 115596 115326
+          115055 114784 114514 114243 113972 113700 113429 113158 112887 112615 112344 112073
+          111802 111530 111259 110988 110716 110445 110174 109903 109631 109360 109089 108817
+          108546 108275 108003 107732 107461 107190 106918 106647 106376 106104 105833 105562
+          105290 105019 104748 104476 104205 103934 103663 103391 103120 102849 102577 102306
+          102035 101763 101492 101221 100950 100678 100407 100136
+        </DataArray>
+        <DataArray type="Float32" Name="pc" 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
+        </DataArray>
+        <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
+          998.631 998.631 998.631 998.631 998.631 998.631 998.632 998.634 998.637 998.642 998.649 998.658
+          998.668 998.678 998.69 998.702 998.715 998.727 998.738 998.749 998.758 998.766 998.773 998.779
+          998.783 998.786 998.789 998.791 998.792 998.793 998.794 998.794 998.794 998.794 998.794 998.794
+          998.794 998.794 998.794 998.794 998.794 998.794 998.794 998.794 998.794 998.793 998.793 998.793
+          998.793 998.793 998.793 998.793 998.793 998.792 998.792 998.792 998.792 998.792 998.792 998.792
+          998.792 998.791 998.791 998.791 998.791 998.791 998.791 998.791 998.791 998.79 998.79 998.79
+          998.79 998.79 998.79 998.79 998.79 998.789 998.789 998.789
+        </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
+        </DataArray>
+        <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
+          946.105 946.106 946.105 946.098 946.075 946.018 945.9 945.683 945.321 944.771 943.999 942.988
+          941.668 939.993 938.14 936.18 934.194 932.256 930.433 928.771 927.303 926.041 924.985 924.123
+          923.435 922.899 922.488 922.181 921.954 921.79 921.674 921.592 921.536 921.498 921.472 921.455
+          921.444 921.437 921.432 921.429 921.427 921.426 921.425 921.425 921.424 921.424 921.424 921.424
+          921.423 921.423 921.423 921.423 921.423 921.423 921.423 921.422 921.422 921.422 921.422 921.422
+          921.422 921.422 921.421 921.421 921.421 921.421 921.421 921.421 921.42 921.42 921.42 921.42
+          921.42 921.42 921.42 921.419 921.419 921.419 921.419 921.418
+        </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
+        </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
+          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
+          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
+          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">
+          291 291 291 291 290.999 290.996 290.991 290.981 290.965 290.941 290.906 290.861
+          290.805 290.74 290.667 290.591 290.512 290.436 290.363 290.297 290.238 290.188 290.145 290.11
+          290.083 290.061 290.044 290.032 290.023 290.016 290.011 290.008 290.006 290.004 290.003 290.002
+          290.002 290.002 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001
+          290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001
+          290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001
+          290.001 290.001 290.001 290.001 290.001 290.001 290.001 290.001
+        </DataArray>
+        <DataArray type="Float32" Name="velocity" NumberOfComponents="3" format="ascii">
+          5e-05 0 0 0.0001 0 0 0.0001 0 0 0.0001 0 0
+          0.0001 0 0 0.0001 0 0 0.0001 0 0 0.0001 0 0
+          9.99999e-05 0 0 9.99997e-05 0 0 9.99995e-05 0 0 9.99991e-05 0 0
+          9.99987e-05 0 0 9.99982e-05 0 0 9.99977e-05 0 0 9.99972e-05 0 0
+          9.99967e-05 0 0 9.99962e-05 0 0 9.99956e-05 0 0 9.99952e-05 0 0
+          9.99948e-05 0 0 9.99944e-05 0 0 9.99941e-05 0 0 9.99938e-05 0 0
+          9.99936e-05 0 0 9.99934e-05 0 0 9.99933e-05 0 0 9.99932e-05 0 0
+          9.99932e-05 0 0 9.99931e-05 0 0 9.99931e-05 0 0 9.99931e-05 0 0
+          9.99931e-05 0 0 9.99931e-05 0 0 9.99931e-05 0 0 9.99931e-05 0 0
+          9.99931e-05 0 0 9.99931e-05 0 0 9.99931e-05 0 0 9.99931e-05 0 0
+          9.99932e-05 0 0 9.99932e-05 0 0 9.99932e-05 0 0 9.99932e-05 0 0
+          9.99932e-05 0 0 9.99932e-05 0 0 9.99932e-05 0 0 9.99933e-05 0 0
+          9.99933e-05 0 0 9.99933e-05 0 0 9.99933e-05 0 0 9.99933e-05 0 0
+          9.99933e-05 0 0 9.99933e-05 0 0 9.99933e-05 0 0 9.99934e-05 0 0
+          9.99934e-05 0 0 9.99934e-05 0 0 9.99934e-05 0 0 9.99934e-05 0 0
+          9.99934e-05 0 0 9.99934e-05 0 0 9.99934e-05 0 0 9.99935e-05 0 0
+          9.99935e-05 0 0 9.99935e-05 0 0 9.99935e-05 0 0 9.99935e-05 0 0
+          9.99935e-05 0 0 9.99935e-05 0 0 9.99936e-05 0 0 9.99936e-05 0 0
+          9.99936e-05 0 0 9.99936e-05 0 0 9.99936e-05 0 0 9.99936e-05 0 0
+          9.99936e-05 0 0 9.99936e-05 0 0 9.99937e-05 0 0 9.99937e-05 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="process rank" 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
+        </DataArray>
+        <DataArray type="Float32" Name="temperatureExact" NumberOfComponents="1" format="ascii">
+          291 291 291 291 291 291 291 291 291 291 291 291
+          291 291 291 291 291 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290 290 290 290 290
+          290 290 290 290 290 290 290 290
+        </DataArray>
+      </CellData>
+      <Points>
+        <DataArray type="Float32" Name="Coordinates" NumberOfComponents="3" format="ascii">
+          0 0 0 0.25 0 0 0 1 0 0.25 1 0
+          0.5 0 0 0.5 1 0 0.75 0 0 0.75 1 0
+          1 0 0 1 1 0 1.25 0 0 1.25 1 0
+          1.5 0 0 1.5 1 0 1.75 0 0 1.75 1 0
+          2 0 0 2 1 0 2.25 0 0 2.25 1 0
+          2.5 0 0 2.5 1 0 2.75 0 0 2.75 1 0
+          3 0 0 3 1 0 3.25 0 0 3.25 1 0
+          3.5 0 0 3.5 1 0 3.75 0 0 3.75 1 0
+          4 0 0 4 1 0 4.25 0 0 4.25 1 0
+          4.5 0 0 4.5 1 0 4.75 0 0 4.75 1 0
+          5 0 0 5 1 0 5.25 0 0 5.25 1 0
+          5.5 0 0 5.5 1 0 5.75 0 0 5.75 1 0
+          6 0 0 6 1 0 6.25 0 0 6.25 1 0
+          6.5 0 0 6.5 1 0 6.75 0 0 6.75 1 0
+          7 0 0 7 1 0 7.25 0 0 7.25 1 0
+          7.5 0 0 7.5 1 0 7.75 0 0 7.75 1 0
+          8 0 0 8 1 0 8.25 0 0 8.25 1 0
+          8.5 0 0 8.5 1 0 8.75 0 0 8.75 1 0
+          9 0 0 9 1 0 9.25 0 0 9.25 1 0
+          9.5 0 0 9.5 1 0 9.75 0 0 9.75 1 0
+          10 0 0 10 1 0 10.25 0 0 10.25 1 0
+          10.5 0 0 10.5 1 0 10.75 0 0 10.75 1 0
+          11 0 0 11 1 0 11.25 0 0 11.25 1 0
+          11.5 0 0 11.5 1 0 11.75 0 0 11.75 1 0
+          12 0 0 12 1 0 12.25 0 0 12.25 1 0
+          12.5 0 0 12.5 1 0 12.75 0 0 12.75 1 0
+          13 0 0 13 1 0 13.25 0 0 13.25 1 0
+          13.5 0 0 13.5 1 0 13.75 0 0 13.75 1 0
+          14 0 0 14 1 0 14.25 0 0 14.25 1 0
+          14.5 0 0 14.5 1 0 14.75 0 0 14.75 1 0
+          15 0 0 15 1 0 15.25 0 0 15.25 1 0
+          15.5 0 0 15.5 1 0 15.75 0 0 15.75 1 0
+          16 0 0 16 1 0 16.25 0 0 16.25 1 0
+          16.5 0 0 16.5 1 0 16.75 0 0 16.75 1 0
+          17 0 0 17 1 0 17.25 0 0 17.25 1 0
+          17.5 0 0 17.5 1 0 17.75 0 0 17.75 1 0
+          18 0 0 18 1 0 18.25 0 0 18.25 1 0
+          18.5 0 0 18.5 1 0 18.75 0 0 18.75 1 0
+          19 0 0 19 1 0 19.25 0 0 19.25 1 0
+          19.5 0 0 19.5 1 0 19.75 0 0 19.75 1 0
+          20 0 0 20 1 0
+        </DataArray>
+      </Points>
+      <Cells>
+        <DataArray type="Int32" Name="connectivity" NumberOfComponents="1" format="ascii">
+          0 1 3 2 1 4 5 3 4 6 7 5
+          6 8 9 7 8 10 11 9 10 12 13 11
+          12 14 15 13 14 16 17 15 16 18 19 17
+          18 20 21 19 20 22 23 21 22 24 25 23
+          24 26 27 25 26 28 29 27 28 30 31 29
+          30 32 33 31 32 34 35 33 34 36 37 35
+          36 38 39 37 38 40 41 39 40 42 43 41
+          42 44 45 43 44 46 47 45 46 48 49 47
+          48 50 51 49 50 52 53 51 52 54 55 53
+          54 56 57 55 56 58 59 57 58 60 61 59
+          60 62 63 61 62 64 65 63 64 66 67 65
+          66 68 69 67 68 70 71 69 70 72 73 71
+          72 74 75 73 74 76 77 75 76 78 79 77
+          78 80 81 79 80 82 83 81 82 84 85 83
+          84 86 87 85 86 88 89 87 88 90 91 89
+          90 92 93 91 92 94 95 93 94 96 97 95
+          96 98 99 97 98 100 101 99 100 102 103 101
+          102 104 105 103 104 106 107 105 106 108 109 107
+          108 110 111 109 110 112 113 111 112 114 115 113
+          114 116 117 115 116 118 119 117 118 120 121 119
+          120 122 123 121 122 124 125 123 124 126 127 125
+          126 128 129 127 128 130 131 129 130 132 133 131
+          132 134 135 133 134 136 137 135 136 138 139 137
+          138 140 141 139 140 142 143 141 142 144 145 143
+          144 146 147 145 146 148 149 147 148 150 151 149
+          150 152 153 151 152 154 155 153 154 156 157 155
+          156 158 159 157 158 160 161 159
+        </DataArray>
+        <DataArray type="Int32" Name="offsets" NumberOfComponents="1" format="ascii">
+          4 8 12 16 20 24 28 32 36 40 44 48
+          52 56 60 64 68 72 76 80 84 88 92 96
+          100 104 108 112 116 120 124 128 132 136 140 144
+          148 152 156 160 164 168 172 176 180 184 188 192
+          196 200 204 208 212 216 220 224 228 232 236 240
+          244 248 252 256 260 264 268 272 276 280 284 288
+          292 296 300 304 308 312 316 320
+        </DataArray>
+        <DataArray type="UInt8" Name="types" NumberOfComponents="1" format="ascii">
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9 9 9 9 9
+          9 9 9 9 9 9 9 9
+        </DataArray>
+      </Cells>
+    </Piece>
+  </UnstructuredGrid>
+</VTKFile>