diff --git a/applications/solvers/incompressible/RichardsFoam/Make/files b/applications/solvers/incompressible/RichardsFoam/Make/files
new file mode 100755
index 000000000..f2b8d52a1
--- /dev/null
+++ b/applications/solvers/incompressible/RichardsFoam/Make/files
@@ -0,0 +1,3 @@
+RichardsFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/RichardsFoam
diff --git a/applications/solvers/incompressible/RichardsFoam/Make/options b/applications/solvers/incompressible/RichardsFoam/Make/options
new file mode 100755
index 000000000..fa15f1245
--- /dev/null
+++ b/applications/solvers/incompressible/RichardsFoam/Make/options
@@ -0,0 +1,5 @@
+EXE_INC = \
+ -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+ -lfiniteVolume
diff --git a/applications/solvers/incompressible/RichardsFoam/RichardsFoam.C b/applications/solvers/incompressible/RichardsFoam/RichardsFoam.C
new file mode 100644
index 000000000..95b45b72b
--- /dev/null
+++ b/applications/solvers/incompressible/RichardsFoam/RichardsFoam.C
@@ -0,0 +1,145 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM 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 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM 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 OpenFOAM. If not, see .
+
+Application
+ RichardsFoam
+
+Description
+ Transient solver for flow in unsaturated porous media
+ With chord slope formulation of the Richards equation.
+ van Genuchten laws for unsaturated hydraulic properties parametrisation
+ Global computation of the convergence criterium
+ Adaptative time stepping with a stabilisation procedure
+ NB 1: use backward scheme for time discretisation
+ NB 2: use only mesh with constant cell volumes
+
+References
+ version 0.0 (develloped with OpenFOAM 2.0.1)
+ Details may be found in:
+ Orgogozo, L., Renon, N., Soulaine, C., Hénon, F., Tomer, S.K., Labat, D.,
+ Pokrovsky, O.S., Sekhar, M., Ababou, R., Quintard, M., Submitted.
+ Mechanistic modelling of water fluxes at the watershed scale: An open source
+ massively parallel solver for Richards equation.
+ Submitted to Computer Physics Communications.
+
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "pimpleControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+# include "setRootCase.H"
+# include "createTime.H"
+# include "createMesh.H"
+
+// pimpleControl pimple(mesh);
+
+# include "readPicardControls.H"
+# include "createFields.H"
+# include "initContinuityErrs.H"
+# include "createTimeControls.H"
+
+ Info<< "\nStarting time loop\n" << endl;
+
+ // starting of the time loop.
+ while (runTime.loop())
+ {
+ // time step control operations.
+# include "readTimeControls.H"
+# include "setDeltaT.H"
+
+// runTime++;
+
+ Info<< "Time = " << runTime.timeName() << nl << endl;
+
+ // Beginning of the stabilisation loop for the stabilised adaptive time
+ // step procedure.
+ for (int cyc = 0; cyc < nMaxCycle; cyc++)
+ {
+ // Beginning of the Picard loop.
+ for (int pic = 0; pic < nIterPicard; pic++)
+ {
+# include "psiEqn.H"
+ }
+
+ // Exit test for the loop associated with the stabilisation cycles
+ // for the adaptive time step procedure.
+ if (crit < precPicard)
+ {
+ break;
+ }
+ else
+ {
+ Info << "Criterion not reached, restart time loop iteration"
+ << "with a smaller time step / Error = " << crit
+ << nl << endl;
+
+ runTime.setDeltaT((1/tFact)*runTime.deltaTValue());
+
+ Info<< "deltaT = " << runTime.deltaTValue() << endl;
+ }
+ // End of the stabilisation cycles loop.
+ }
+
+ // Warning test in case of convergence failure of the Picard loop.
+ if (crit >= precPicard)
+ {
+ Info<< "Convergence failure / Error = " << crit << nl << endl;
+ currentPicard = nIterPicard;
+ }
+
+ // Final updating of the result fields before going to the next time
+ // iteration.
+ psi_tmp = psi;
+
+ thtil_tmp = 0.5*
+ (
+ (1 + sign(psi_tmp)) + (1 - sign(psi_tmp))*
+ pow((1 + pow(mag(alpha*psi_tmp),n)), - (1 - (1/n)))
+ );
+
+ theta = (thetas - thetar)*thtil + thetar;
+
+ U = - Krel*((fvc::grad(psi)) + vuz);
+
+ // Writting of the result.
+ runTime.write();
+
+ Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+ << " ClockTime = " << runTime.elapsedClockTime() << " s"
+ << nl << endl;
+
+ // end of the time loop.
+ }
+
+ Info<< "End\n" << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/incompressible/RichardsFoam/createFields.H b/applications/solvers/incompressible/RichardsFoam/createFields.H
new file mode 100644
index 000000000..246cc0476
--- /dev/null
+++ b/applications/solvers/incompressible/RichardsFoam/createFields.H
@@ -0,0 +1,209 @@
+Info<< "Reading transportProperties\n" << endl;
+
+// reading of the physical constant associated with the considered problem.
+// Localisation of the data within the considered case:
+// constant/transportProperties
+IOdictionary transportProperties
+(
+ IOobject
+ (
+ "transportProperties",
+ runTime.constant(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::NO_WRITE
+ )
+);
+
+dimensionedScalar K
+(
+ transportProperties.lookup("K")
+);
+
+dimensionedScalar alpha
+(
+ transportProperties.lookup("alpha")
+);
+
+dimensionedScalar thetas
+(
+ transportProperties.lookup("thetas")
+);
+
+dimensionedScalar thetar
+(
+ transportProperties.lookup("thetar")
+);
+
+dimensionedScalar n
+(
+ transportProperties.lookup("n")
+);
+
+dimensionedScalar C
+(
+ transportProperties.lookup("C")
+);
+
+dimensionedScalar S
+(
+ transportProperties.lookup("S")
+);
+
+// declaration of the variable and results fields
+
+// Water velocity field [m/s]
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+ IOobject
+ (
+ "U",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+// Water saturation field [-]
+Info<< "Reading field theta\n" << endl;
+volScalarField theta
+(
+ IOobject
+ (
+ "theta",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+// Water pressure field [m] - field of resolution.
+Info<< "Reading field psi\n" << endl;
+volScalarField psi
+(
+ IOobject
+ (
+ "psi",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+// Field of residuals for the Picard loop [m].
+Info<< "Reading field err\n" << endl;
+volScalarField err
+(
+ IOobject
+ (
+ "err",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+// Dimensionned unit scalar field [m].
+Info<< "Reading field vuz\n" << endl;
+volVectorField vuz
+(
+ IOobject
+ (
+ "vuz",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+// Dimensionless unit vertical upward vector field.
+Info<< "Reading field usf\n" << endl;
+volScalarField usf
+(
+ IOobject
+ (
+ "usf",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+);
+
+# include "createPhi.H"
+
+// Initialisation of the scalar containing the number of mesh cells. Note the
+// use of gSum instead of sum.
+double nbMesh;
+nbMesh = gSum(usf);
+
+// Initialisation of the scalar containing the residual for the exit test of the
+// Picard loop.
+double crit;
+crit=0.;
+
+// Initialisation of the token which counts the number of Picard iteraion for
+// the adaptive time step procedure.
+int currentPicard;
+currentPicard = nIterPicard-3;
+
+// Initialisation of the token which counts the number of Stabilisation cycle
+// for the stabilisation of the adaptive time step procedure.
+int sc;
+sc = 0;
+
+// Initialisation of the field of altitudes.
+volVectorField positionVector = mesh.C();
+volScalarField z = positionVector.component(vector::Z);
+
+// Initialisation of the intermediate fields for the Picard loop.
+volScalarField psi_tmp = psi;
+volScalarField psim1 = psi;
+
+// Initialisation of the varying transport properties for the Picard loop.
+volScalarField thtil =
+ 0.5*
+ (
+ (1 + sign(psi)) + (1 - sign(psi))*
+ pow((1 + pow(mag(alpha*psi),n)), - (1 - (1/n)))
+ );
+
+volScalarField thtil_tmp =
+ 0.5*
+ (
+ (1 + sign(psi_tmp)) + (1-sign(psi_tmp))*
+ pow((1 + pow(mag(alpha*psi_tmp),n)), - (1 - (1/n)))
+ );
+
+volScalarField Krel =
+ 0.5*
+ (
+ (1 + sign(psi))*K + (1 - sign(psi))*K*pow(thtil,0.5)*
+ pow((1 - pow((1 - pow(thtil,(n/(n - 1)))),(1 - (1/n)))),2)
+ );
+
+volScalarField Crel =
+ S + 0.5*
+ (
+ (1 - sign(psi))*((thetas - thetar)*(thtil - thtil_tmp)*
+ (1./((usf*pos(psi - psi_tmp)*pos(psi_tmp - psi)) + psi - psi_tmp)))
+ );
+
+// Initialisation of the gravity term.
+volVectorField gradk = fvc::grad(Krel);
+volScalarField gradkz = gradk.component(vector::Z);
+
+// Initialisation of the velocity field.
+U = - Krel*((fvc::grad(psi)) + vuz);
diff --git a/applications/solvers/incompressible/RichardsFoam/psiEqn.H b/applications/solvers/incompressible/RichardsFoam/psiEqn.H
new file mode 100644
index 000000000..a7f03a9d1
--- /dev/null
+++ b/applications/solvers/incompressible/RichardsFoam/psiEqn.H
@@ -0,0 +1,60 @@
+psim1 = psi;
+psi = psi_tmp;
+
+// Resolution of the linear system.
+fvScalarMatrix psiEqn
+(
+ Crel*fvm::ddt(psi)
+ ==
+ fvm::laplacian(Krel, psi, "laplacian(Krel,psi)")
+ + gradkz
+);
+psiEqn.relax();
+psiEqn.solve();
+
+// Update of the varying transport properties.
+thtil = 0.5*
+(
+ (1 + sign(psi)) + (1 - sign(psi))*
+ pow((1+pow(mag(alpha*psi),n)),-(1-(1/n)))
+);
+Krel = 0.5*
+(
+ (1 + sign(psi))*K + (1 - sign(psi))*K*pow(thtil, 0.5)*
+ pow((1 - pow((1 - pow(thtil, (n/(n - 1)))),(1 - (1/n)))), 2)
+);
+Crel= S + 0.5*
+(
+ (1 - sign(psi))*
+ (
+ (thetas - thetar)*(thtil - thtil_tmp)*
+ (
+ 1./((usf*pos(psi - psi_tmp)*pos(psi_tmp - psi))
+ + psi - psi_tmp)
+ )
+ )
+);
+
+// Update of the gravity term.
+gradk = fvc::grad(Krel);
+gradkz = gradk.component(2);
+
+// Computation of the field of residuals for the exit test of
+// the Picard loop.
+err = psi - psim1;
+
+// Computation of the residual for the exit test of the Picard
+// loop. Note the use of gSum instead of sum.
+crit = gSumMag(err)/nbMesh;
+
+// exit test for the Picard loop.
+if (crit < precPicard)
+{
+ Info << " Erreur = " << crit
+ << " Picard = " << pic
+ << nl << endl;
+
+ currentPicard=pic;
+
+ break;
+}
diff --git a/applications/solvers/incompressible/RichardsFoam/readPicardControls.H b/applications/solvers/incompressible/RichardsFoam/readPicardControls.H
new file mode 100644
index 000000000..68735fad5
--- /dev/null
+++ b/applications/solvers/incompressible/RichardsFoam/readPicardControls.H
@@ -0,0 +1,25 @@
+// Reading of the controling parameters of the Picard loop and of the
+// Stabilisation cycles loop for the adaptive time step procedure.
+// Localisation of the data within the considered case: system/fvSolution
+const dictionary& picardDict = mesh.solutionDict().subDict("Picard");
+
+// Maximum number of Picard iterations.
+const int nIterPicard =
+ picardDict.lookupOrDefault("nIterPicard", 1);
+
+// Maximum number of stabilisation cycles.
+const int nMaxCycle =
+ picardDict.lookupOrDefault("nMaxCycle", 1);
+
+// Number of time iterations with low number of Picard iteration beyond which
+// the time step is increased.
+const int stabilisationThreshold =
+ picardDict.lookupOrDefault("stabilisationThreshold", 1);
+
+// Exit criterion for the Picard loop.
+const double precPicard =
+ picardDict.lookupOrDefault("precPicard", 1.);
+
+// Factor of increase/decrease of the time step.
+const double tFact =
+ picardDict.lookupOrDefault("tFact", 1.);
diff --git a/applications/solvers/incompressible/RichardsFoam/setDeltaT.H b/applications/solvers/incompressible/RichardsFoam/setDeltaT.H
new file mode 100644
index 000000000..c597cda93
--- /dev/null
+++ b/applications/solvers/incompressible/RichardsFoam/setDeltaT.H
@@ -0,0 +1,82 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2004-2010 OpenCFD Ltd.
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM 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 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM 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 OpenFOAM. If not, see .
+
+Global
+ setDeltaT
+
+Description
+ Reset the timestep to maintain a low number of Picard iteration and a good
+ convergence. Reduction of time-step is immediate, but increase is damped to
+ avoid unstable oscillations.
+
+\*---------------------------------------------------------------------------*/
+
+scalar deltaTFact = 1.;
+
+if (adjustTimeStep)
+{
+ deltaTFact = 1.;
+
+ // If There is a too high number of Picard iteration, decrease of the time
+ // step.
+ if (currentPicard > nIterPicard-2)
+ {
+ deltaTFact = 1./tFact;
+ }
+
+ if (currentPicard >= 3)
+ {
+ sc = 0;
+ }
+
+ if ((currentPicard < 3) && (sc < stabilisationThreshold))
+ {
+ sc = sc+1;
+ }
+
+ // If there is more than 'stabilisationThreshold' iteration that the number
+ // of Picard iteration is low, increase of the time step.
+ if ((currentPicard < 3) && (sc == stabilisationThreshold))
+ {
+ deltaTFact = tFact;
+ sc = 0;
+ }
+
+ // Reset of the time step if needed (time step is always lower than
+ // 'maxDeltaT').
+ if (deltaTFact != 1.)
+ {
+ runTime.setDeltaT
+ (
+ min
+ (
+ deltaTFact*runTime.deltaTValue(),
+ maxDeltaT
+ )
+ );
+ }
+
+ Info<< "deltaT = " << runTime.deltaTValue() << endl;
+}
+
+// ************************************************************************* //
diff --git a/applications/solvers/incompressible/RichardsFoam/spatialMeanValueRichardsonFoam/Make/files b/applications/solvers/incompressible/RichardsFoam/spatialMeanValueRichardsonFoam/Make/files
new file mode 100755
index 000000000..32c9117fd
--- /dev/null
+++ b/applications/solvers/incompressible/RichardsFoam/spatialMeanValueRichardsonFoam/Make/files
@@ -0,0 +1,3 @@
+spatialMeanValueRichardsonFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/spatialMeanValueRichardsonFoam
diff --git a/applications/solvers/incompressible/RichardsFoam/spatialMeanValueRichardsonFoam/Make/options b/applications/solvers/incompressible/RichardsFoam/spatialMeanValueRichardsonFoam/Make/options
new file mode 100755
index 000000000..fa15f1245
--- /dev/null
+++ b/applications/solvers/incompressible/RichardsFoam/spatialMeanValueRichardsonFoam/Make/options
@@ -0,0 +1,5 @@
+EXE_INC = \
+ -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+ -lfiniteVolume
diff --git a/applications/solvers/incompressible/RichardsFoam/spatialMeanValueRichardsonFoam/spatialMeanValueRichardsonFoam.C b/applications/solvers/incompressible/RichardsFoam/spatialMeanValueRichardsonFoam/spatialMeanValueRichardsonFoam.C
new file mode 100644
index 000000000..298304b98
--- /dev/null
+++ b/applications/solvers/incompressible/RichardsFoam/spatialMeanValueRichardsonFoam/spatialMeanValueRichardsonFoam.C
@@ -0,0 +1,110 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM 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 3 of the License, or
+ (at your option) any later version.
+
+ OpenFOAM 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 OpenFOAM. If not, see .
+
+Application
+ spatialMeanValue
+
+Description
+ Calculates the spatial mean of the specified volScalarField over the
+ whole domain (for regular grid).
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "cyclicPolyPatch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+// Main program:
+
+int main(int argc, char *argv[])
+{
+# include "addRegionOption.H"
+
+ timeSelector::addOptions();
+ argList::validArgs.append("fieldName");
+
+# include "setRootCase.H"
+# include "createTime.H"
+
+ instantList timeDirs = timeSelector::select0(runTime, args);
+
+# include "createNamedMesh.H"
+
+ word fieldName(args.additionalArgs()[0]);
+
+ Info << "Calculating avarage pressure heads:" << endl;
+
+ forAll(timeDirs, timeI)
+ {
+ runTime.setTime(timeDirs[timeI], timeI);
+
+ IOobject fieldHeader
+ (
+ fieldName,
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ
+ );
+
+ // Check field exists
+ if (fieldHeader.headerOk())
+ {
+ mesh.readUpdate();
+
+ // Read field and calc mean, for regular grid
+ if (fieldHeader.headerClassName() == volScalarField::typeName)
+ {
+ volScalarField field(fieldHeader, mesh);
+
+ int nbMesh;
+ nbMesh = 0;
+
+ forAll(field, cellI)
+ {
+ nbMesh++;
+ }
+
+ Info<< runTime.timeName()<< " "
+ << sum(field).value()/nbMesh<< " "
+ << endl;
+ }
+ else
+ {
+ FatalError
+ << "Only possible for volScalarField "<< "s "
+ << nl << exit(FatalError);
+ }
+ }
+ else
+ {
+ Info<< " No field " << fieldName << endl;
+ }
+ }
+
+ Info<< "End\n" << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/U b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/U
new file mode 100755
index 000000000..48e65a924
--- /dev/null
+++ b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/U
@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: 1.7.x |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class volVectorField;
+ location "0";
+ object U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 1 -1 0 0 0 0];
+
+internalField uniform (0 0 0);
+
+boundaryField
+{
+ top
+ {
+ type zeroGradient;
+ }
+ bottom
+ {
+ type zeroGradient;
+ }
+ sides
+ {
+ type empty;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/err b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/err
new file mode 100755
index 000000000..9dc130b23
--- /dev/null
+++ b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/err
@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: 1.7.x |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class volScalarField;
+ location "0";
+ object err;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 1 0 0 0 0 0];
+
+internalField uniform 1;
+
+boundaryField
+{
+ top
+ {
+ type zeroGradient;
+ }
+ bottom
+ {
+ type zeroGradient;
+ }
+ sides
+ {
+ type empty;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/psi b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/psi
new file mode 100755
index 000000000..13b5ae646
--- /dev/null
+++ b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/psi
@@ -0,0 +1,43 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: 1.7.x |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class volScalarField;
+ location "0";
+ object psi;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 1 0 0 0 0 0];
+
+internalField uniform -1;
+
+boundaryField
+{
+ top
+ {
+ type fixedValue;
+ value uniform 0.01;
+ }
+
+ bottom
+ {
+ type fixedGradient;
+ gradient uniform 0;
+ }
+
+ sides
+ {
+ type empty;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/theta b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/theta
new file mode 100755
index 000000000..ce23341e4
--- /dev/null
+++ b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/theta
@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: 1.7.x |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class volScalarField;
+ location "0";
+ object theta;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 0 0 0 0 0 0];
+
+internalField uniform 0;
+
+boundaryField
+{
+ top
+ {
+ type zeroGradient;
+ }
+ bottom
+ {
+ type zeroGradient;
+ }
+ sides
+ {
+ type empty;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/usf b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/usf
new file mode 100755
index 000000000..fd1f12cb9
--- /dev/null
+++ b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/usf
@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: 1.7.x |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class volScalarField;
+ location "0";
+ object usf;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 1 0 0 0 0 0];
+
+internalField uniform 1.;
+
+boundaryField
+{
+ top
+ {
+ type zeroGradient;
+ }
+ bottom
+ {
+ type zeroGradient;
+ }
+ sides
+ {
+ type empty;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/vuz b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/vuz
new file mode 100755
index 000000000..6601b9577
--- /dev/null
+++ b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/0/vuz
@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: 1.7.x |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class volVectorField;
+ location "0";
+ object vuz;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions [0 0 0 0 0 0 0];
+
+internalField uniform (0 0 1);
+
+boundaryField
+{
+ top
+ {
+ type zeroGradient;
+ }
+ bottom
+ {
+ type zeroGradient;
+ }
+ sides
+ {
+ type empty;
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/Allclean b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/Allclean
new file mode 100755
index 000000000..e80d74cac
--- /dev/null
+++ b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/Allclean
@@ -0,0 +1,7 @@
+#!/bin/sh
+
+# Source tutorial clean functions
+. $WM_PROJECT_DIR/bin/tools/CleanFunctions
+
+cleanCase
+
diff --git a/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/Allrun b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/Allrun
new file mode 100755
index 000000000..cfb31d809
--- /dev/null
+++ b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/Allrun
@@ -0,0 +1,7 @@
+#!/bin/sh
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+runApplication blockMesh
+runApplication RichardsFoam
+runApplication spatialMeanValueRichardsonFoam psi
diff --git a/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/constant/polyMesh/blockMeshDict b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/constant/polyMesh/blockMeshDict
new file mode 100644
index 000000000..61543cfb7
--- /dev/null
+++ b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/constant/polyMesh/blockMeshDict
@@ -0,0 +1,68 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: 1.7.1 |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ object blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1;
+
+vertices
+(
+ (0 0 0)
+ (1 0 0)
+ (1 1 0)
+ (0 1 0)
+
+ (0 0 1)
+ (1 0 1)
+ (1 1 1)
+ (0 1 1)
+);
+
+blocks
+(
+ hex (0 1 2 3 4 5 6 7) (1 1 100) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+patches
+
+(
+
+ patch top
+ (
+ (4 5 6 7)
+ )
+
+ patch bottom
+ (
+ (0 1 2 3)
+ )
+
+ empty sides
+ (
+ (1 2 6 5)
+ (3 2 6 7)
+ (0 3 7 4)
+ (0 1 5 4)
+ )
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/constant/polyMesh/boundary b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/constant/polyMesh/boundary
new file mode 100644
index 000000000..340034118
--- /dev/null
+++ b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/constant/polyMesh/boundary
@@ -0,0 +1,41 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | foam-extend: Open Source CFD |
+| \\ / O peration | Version: 4.0 |
+| \\ / A nd | Web: http://www.foam-extend.org |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class polyBoundaryMesh;
+ location "constant/polyMesh";
+ object boundary;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+3
+(
+ top
+ {
+ type patch;
+ nFaces 1;
+ startFace 99;
+ }
+ bottom
+ {
+ type patch;
+ nFaces 1;
+ startFace 100;
+ }
+ sides
+ {
+ type empty;
+ nFaces 400;
+ startFace 101;
+ }
+)
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/constant/transportProperties b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/constant/transportProperties
new file mode 100755
index 000000000..ddeca6530
--- /dev/null
+++ b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/constant/transportProperties
@@ -0,0 +1,32 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: 1.7.1 |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ location "constant";
+ object transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+
+
+K K [ 0 1 -1 0 0 0 0 ] 0.00000288889;
+alpha alpha [ 0 -1 0 0 0 0 0 ] 3.6;
+n n [ 0 0 0 0 0 0 0 ] 1.56;
+thetas thetas [ 0 0 0 0 0 0 0 ] 0.43;
+thetar thetar [ 0 0 0 0 0 0 0 ] 0.078;
+S S [ 0 -1 0 0 0 0 0 ] 0.00001;
+C C [ 0 -1 0 0 0 0 0 ] 1;
+
+
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/system/controlDict b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/system/controlDict
new file mode 100644
index 000000000..d5393573a
--- /dev/null
+++ b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/system/controlDict
@@ -0,0 +1,78 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: 1.7.1 |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ location "system";
+ object controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application porousUnsaturatedFoam;
+
+startFrom startTime;
+
+startTime 0;
+
+stopAt endTime;
+
+endTime 86400;
+
+deltaT 300;
+
+writeControl adjustableRunTime;
+
+writeInterval 7200;
+
+purgeWrite 0;
+
+writeFormat ascii;
+
+writePrecision 6;
+
+writeCompression uncompressed;
+
+timeFormat general;
+
+timePrecision 6;
+
+runTimeModifiable yes;
+
+adjustTimeStep yes;
+
+maxDeltaT 3600;
+
+functions
+ (
+ probes1
+ {
+ type probes; // Type of functionObject
+ // Where to load it from (if not already in solver)
+ functionObjectLibs ("libsampling.so");
+ probeLocations // Locations to be probed. runTime modifiable!
+ (
+ (0.5 0.5 1.)
+ );
+ // Fields to be probed. runTime modifiable!
+ fields
+ (
+ U
+ );
+
+ outputControl timeStep;
+ outputInterval 1;
+ }
+ );
+
+
+//libs ( "libOpenFOAM.so" );
+//libs ("libuserBCs.so");
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/system/fvSchemes b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/system/fvSchemes
new file mode 100755
index 000000000..74173cf7e
--- /dev/null
+++ b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/system/fvSchemes
@@ -0,0 +1,58 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: 1.7.1 |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ location "system";
+ object fvSchemes;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+ddtSchemes
+{
+ default backward;
+}
+
+gradSchemes
+{
+ default Gauss linear;
+
+}
+
+divSchemes
+{
+ default none;
+}
+
+laplacianSchemes
+{
+ default none;
+
+laplacian(Krel,psi) Gauss linear corrected;
+}
+
+interpolationSchemes
+{
+ default Gauss vanLeer;
+}
+
+snGradSchemes
+{
+ default Gauss vanLeer;
+}
+
+fluxRequired
+{
+ default yes;
+ U;
+}
+
+
+// ************************************************************************* //
diff --git a/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/system/fvSolution b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/system/fvSolution
new file mode 100755
index 000000000..6b4a77193
--- /dev/null
+++ b/tutorials/incompressible/RichardsonFoam/oneDimensionalSoilTransport/system/fvSolution
@@ -0,0 +1,44 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| ========= | |
+| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
+| \\ / O peration | Version: 1.7.1 |
+| \\ / A nd | Web: www.OpenFOAM.com |
+| \\/ M anipulation | |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+ version 2.0;
+ format ascii;
+ class dictionary;
+ location "system";
+ object fvSolution;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+solvers
+{
+
+ psi
+ {
+ solver PCG;
+ preconditioner DIC;
+ tolerance 1e-6;
+ relTol 0;
+ }
+}
+
+Picard
+{
+ nIterPicard 10; //must be strictly higher than 5 ; 10 should be the best value//
+ nMaxCycle 10;
+ stabilisationThreshold 10;
+ precPicard 1e-4;
+ tFact 1.3; //1 means fixed time step ; should be lower than 1.5//
+}
+
+relaxationFactor
+{
+// psi 0.3;
+};
+
+// ************************************************************************* //