diff --git a/Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/c_RelativePermeabilityVanGenuchten.md b/Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/c_RelativePermeabilityVanGenuchten.md
new file mode 100644
index 0000000000000000000000000000000000000000..aa3bffae1ae88769d43e7c20296d2fb4788e02ed
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/c_RelativePermeabilityVanGenuchten.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::RelativePermeabilityVanGenuchten
diff --git a/Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/t_exponent.md b/Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/t_exponent.md
new file mode 100644
index 0000000000000000000000000000000000000000..0f0ec649ba83f9c2a0c3f7c2889980c807d1edda
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/t_exponent.md
@@ -0,0 +1 @@
+The exponent of the van Genuchten relative permeability function.
diff --git a/Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/t_minimum_relative_permeability_liquid.md b/Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/t_minimum_relative_permeability_liquid.md
new file mode 100644
index 0000000000000000000000000000000000000000..926d443546280f308722c4265ec92e2cf8703deb
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/t_minimum_relative_permeability_liquid.md
@@ -0,0 +1 @@
+The minimum relative permeability of the liquid (or wetting) phase.
diff --git a/Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/t_residual_gas_saturation.md b/Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/t_residual_gas_saturation.md
new file mode 100644
index 0000000000000000000000000000000000000000..d8631eb27f4627345443d3cb426330a84c325db1
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/t_residual_gas_saturation.md
@@ -0,0 +1 @@
+The smallest degree of saturation of the gas (or non-wetting) phase in the medium.
diff --git a/Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/t_residual_liquid_saturation.md b/Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/t_residual_liquid_saturation.md
new file mode 100644
index 0000000000000000000000000000000000000000..95ce0ed4b319db3ecc0cfc87b42c611dcd572210
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/RelativePermeabilityVanGenuchten/t_residual_liquid_saturation.md
@@ -0,0 +1 @@
+The smallest degree of saturation of the liquid (or wetting) phase in the medium.
diff --git a/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/c_SaturationVanGenuchten.md b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/c_SaturationVanGenuchten.md
new file mode 100644
index 0000000000000000000000000000000000000000..2a8f10769050f001e0532878f196e583db65c309
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/c_SaturationVanGenuchten.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::SaturationVanGenuchten
diff --git a/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_entry_pressure.md b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_entry_pressure.md
new file mode 100644
index 0000000000000000000000000000000000000000..e3d246498cefcd2402e2df304f0047534317a04a
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_entry_pressure.md
@@ -0,0 +1 @@
+The required pressure for a non-wetting fluid to enter a porous medium (the capillary pressure at full saturation).
diff --git a/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_exponent.md b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_exponent.md
new file mode 100644
index 0000000000000000000000000000000000000000..260b1dd015e7c9cb8b04a65e4df7587dc81b9224
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_exponent.md
@@ -0,0 +1 @@
+The exponent of the van Genuchten saturation function.
diff --git a/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_residual_gas_saturation.md b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_residual_gas_saturation.md
new file mode 100644
index 0000000000000000000000000000000000000000..d8631eb27f4627345443d3cb426330a84c325db1
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_residual_gas_saturation.md
@@ -0,0 +1 @@
+The smallest degree of saturation of the gas (or non-wetting) phase in the medium.
diff --git a/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_residual_liquid_saturation.md b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_residual_liquid_saturation.md
new file mode 100644
index 0000000000000000000000000000000000000000..95ce0ed4b319db3ecc0cfc87b42c611dcd572210
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/SaturationVanGenuchten/t_residual_liquid_saturation.md
@@ -0,0 +1 @@
+The smallest degree of saturation of the liquid (or wetting) phase in the medium.
diff --git a/MaterialLib/MPL/CreateProperty.cpp b/MaterialLib/MPL/CreateProperty.cpp
index abdb35f6433a1175357f7c93b7b55db84aa27b0d..a31c10b166ab1168cef39dce02bf82324657c490 100644
--- a/MaterialLib/MPL/CreateProperty.cpp
+++ b/MaterialLib/MPL/CreateProperty.cpp
@@ -79,6 +79,16 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
return createRelPermLiakopoulos(config);
}
+ if (boost::iequals(property_type, "SaturationVanGenuchten"))
+ {
+ return createSaturationVanGenuchten(config);
+ }
+
+ if (boost::iequals(property_type, "RelativePermeabilityVanGenuchten"))
+ {
+ return createRelPermVanGenuchten(config);
+ }
+
// If none of the above property types are found, OGS throws an error.
OGS_FATAL("The specified component property type '%s' was not recognized",
property_type.c_str());
diff --git a/MaterialLib/MPL/Properties/CreateProperties.h b/MaterialLib/MPL/Properties/CreateProperties.h
index ec8c0508cea1e14266fa167df8a53197ebd159c3..db967698b4f1fbdd8e98cbc979e9110818fef046 100644
--- a/MaterialLib/MPL/Properties/CreateProperties.h
+++ b/MaterialLib/MPL/Properties/CreateProperties.h
@@ -18,5 +18,7 @@
#include "CreateParameterProperty.h"
#include "CreateRelPermBrooksCorey.h"
#include "CreateRelPermLiakopoulos.h"
+#include "CreateRelPermVanGenuchten.h"
#include "CreateSaturationBrooksCorey.h"
-#include "CreateSaturationLiakopoulos.h"
\ No newline at end of file
+#include "CreateSaturationLiakopoulos.h"
+#include "CreateSaturationVanGenuchten.h"
diff --git a/MaterialLib/MPL/Properties/CreateRelPermVanGenuchten.cpp b/MaterialLib/MPL/Properties/CreateRelPermVanGenuchten.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..332f0180abaeda85ff32b179a16746162f505ea8
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreateRelPermVanGenuchten.cpp
@@ -0,0 +1,47 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "BaseLib/ConfigTree.h"
+#include "RelPermVanGenuchten.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<RelPermVanGenuchten> createRelPermVanGenuchten(
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{properties__property__type}
+ config.checkConfigParameter("type", "RelativePermeabilityVanGenuchten");
+ DBUG("Create RelativePermeabilityVanGenuchten medium property");
+
+ auto const residual_liquid_saturation =
+ //! \ogs_file_param{properties__property__RelativePermeabilityVanGenuchten__residual_liquid_saturation}
+ config.getConfigParameter<double>("residual_liquid_saturation");
+ auto const residual_gas_saturation =
+ //! \ogs_file_param{properties__property__RelativePermeabilityVanGenuchten__residual_gas_saturation}
+ config.getConfigParameter<double>("residual_gas_saturation");
+ auto const min_relative_permeability_liquid =
+ //! \ogs_file_param{properties__property__RelativePermeabilityVanGenuchten__minimum_relative_permeability_liquid}
+ config.getConfigParameter<double>(
+ "minimum_relative_permeability_liquid");
+ auto const exponent =
+ //! \ogs_file_param{properties__property__RelativePermeabilityVanGenuchten__exponent}
+ config.getConfigParameter<double>("exponent");
+ if (exponent <= 0. || exponent >= 1.)
+ {
+ OGS_FATAL("Exponent must be in the (0, 1) range.");
+ }
+
+ return std::make_unique<RelPermVanGenuchten>(
+ residual_liquid_saturation,
+ residual_gas_saturation,
+ min_relative_permeability_liquid,
+ exponent);
+}
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CreateRelPermVanGenuchten.h b/MaterialLib/MPL/Properties/CreateRelPermVanGenuchten.h
new file mode 100644
index 0000000000000000000000000000000000000000..046b31552bdbcd42fb4e621b1049e815b02d271d
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreateRelPermVanGenuchten.h
@@ -0,0 +1,29 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+class RelPermVanGenuchten;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<RelPermVanGenuchten> createRelPermVanGenuchten(
+ BaseLib::ConfigTree const& config);
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CreateSaturationVanGenuchten.cpp b/MaterialLib/MPL/Properties/CreateSaturationVanGenuchten.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9b77c416906b13e37db49ca627fa1fcb1c7321e5
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreateSaturationVanGenuchten.cpp
@@ -0,0 +1,41 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "BaseLib/ConfigTree.h"
+#include "SaturationVanGenuchten.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<SaturationVanGenuchten> createSaturationVanGenuchten(
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{properties__property__type}
+ config.checkConfigParameter("type", "SaturationVanGenuchten");
+
+ DBUG("Create SaturationVanGenuchten medium property");
+
+ auto const residual_liquid_saturation =
+ //! \ogs_file_param{properties__property__SaturationVanGenuchten__residual_liquid_saturation}
+ config.getConfigParameter<double>("residual_liquid_saturation");
+ auto const residual_gas_saturation =
+ //! \ogs_file_param{properties__property__SaturationVanGenuchten__residual_gas_saturation}
+ config.getConfigParameter<double>("residual_gas_saturation");
+ auto const exponent =
+ //! \ogs_file_param{properties__property__SaturationVanGenuchten__exponent}
+ config.getConfigParameter<double>("exponent");
+ auto const entry_pressure =
+ //! \ogs_file_param{properties__property__SaturationVanGenuchten__entry_pressure}
+ config.getConfigParameter<double>("entry_pressure");
+
+ return std::make_unique<SaturationVanGenuchten>(residual_liquid_saturation,
+ residual_gas_saturation,
+ exponent, entry_pressure);
+}
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CreateSaturationVanGenuchten.h b/MaterialLib/MPL/Properties/CreateSaturationVanGenuchten.h
new file mode 100644
index 0000000000000000000000000000000000000000..01ffabffe1a515a902b8e95c663136704b72495b
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreateSaturationVanGenuchten.h
@@ -0,0 +1,29 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+class SaturationVanGenuchten;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<SaturationVanGenuchten> createSaturationVanGenuchten(
+ BaseLib::ConfigTree const& config);
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/Properties.h b/MaterialLib/MPL/Properties/Properties.h
index 271a4966fb4c48d575a8b5dc801114c0a9a59892..eab18c0c383aaae33ada9cf9b1c29d4fad3d63da 100644
--- a/MaterialLib/MPL/Properties/Properties.h
+++ b/MaterialLib/MPL/Properties/Properties.h
@@ -18,5 +18,7 @@
#include "ParameterProperty.h"
#include "RelPermBrooksCorey.h"
#include "RelPermLiakopoulos.h"
+#include "RelPermVanGenuchten.h"
#include "SaturationBrooksCorey.h"
#include "SaturationLiakopoulos.h"
+#include "SaturationVanGenuchten.h"
diff --git a/MaterialLib/MPL/Properties/RelPermVanGenuchten.cpp b/MaterialLib/MPL/Properties/RelPermVanGenuchten.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..efc74563608b2203de4f4e79d90531c3acdb6a2b
--- /dev/null
+++ b/MaterialLib/MPL/Properties/RelPermVanGenuchten.cpp
@@ -0,0 +1,95 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "RelPermVanGenuchten.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include "MaterialLib/MPL/Medium.h"
+
+namespace MaterialPropertyLib
+{
+RelPermVanGenuchten::RelPermVanGenuchten(
+ double const residual_liquid_saturation,
+ double const residual_gas_saturation,
+ double const min_relative_permeability_liquid,
+ double const exponent)
+ : _S_L_res(residual_liquid_saturation),
+ _S_L_max(1. - residual_gas_saturation),
+ _k_rel_min(min_relative_permeability_liquid),
+ _m(exponent)
+{
+ if (!(_m > 0 && _m < 1))
+ {
+ OGS_FATAL(
+ "The exponent value m = %g of van Genuchten relative permeability "
+ "model, is out of its range of (0, 1)",
+ _m);
+ }
+}
+
+PropertyDataType RelPermVanGenuchten::value(
+ VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& /*pos*/,
+ double const /*t*/) const
+{
+ double const S_L = std::clamp(
+ std::get<double>(
+ variable_array[static_cast<int>(Variable::liquid_saturation)]),
+ _S_L_res, _S_L_max);
+
+ double const S_eff = (S_L - _S_L_res) / (_S_L_max - _S_L_res);
+ double const v = 1. - std::pow(1. - std::pow(S_eff, 1. / _m), _m);
+ double const k_rel = std::sqrt(S_eff) * v * v;
+ return std::max(_k_rel_min, k_rel);
+}
+
+PropertyDataType RelPermVanGenuchten::dValue(
+ VariableArray const& variable_array, Variable const primary_variable,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/) const
+{
+ (void)primary_variable;
+ assert((primary_variable == Variable::liquid_saturation) &&
+ "RelativePermeabilityVanGenuchten::dValue is implemented for "
+ "derivatives with respect to liquid saturation only.");
+ double const S_L = std::clamp(
+ std::get<double>(
+ variable_array[static_cast<int>(Variable::liquid_saturation)]),
+ _S_L_res, _S_L_max);
+
+ double const S_eff = (S_L - _S_L_res) / (_S_L_max - _S_L_res);
+ if (S_eff <= 0) // prevent division by zero
+ {
+ return 0;
+ }
+
+ if (S_eff >= 1) // prevent taking root of zero
+ {
+ return 0;
+ }
+
+ double const S_eff_to_1_over_m = std::pow(S_eff, 1. / _m);
+ double const v = 1. - std::pow(1. - S_eff_to_1_over_m, _m);
+ double const sqrt_S_eff = std::sqrt(S_eff);
+ double const k_rel = sqrt_S_eff * v * v;
+
+ if (k_rel < _k_rel_min)
+ {
+ return 0;
+ }
+
+ return (0.5 * v * v / sqrt_S_eff +
+ 2. * sqrt_S_eff * v * std::pow(1. - S_eff_to_1_over_m, _m - 1.) *
+ S_eff_to_1_over_m / S_eff) /
+ (_S_L_max - _S_L_res);
+}
+
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/RelPermVanGenuchten.h b/MaterialLib/MPL/Properties/RelPermVanGenuchten.h
new file mode 100644
index 0000000000000000000000000000000000000000..a317527cad8068cc1471bf4569ce694485df1c1c
--- /dev/null
+++ b/MaterialLib/MPL/Properties/RelPermVanGenuchten.h
@@ -0,0 +1,66 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+#pragma once
+
+#include "BaseLib/ConfigTree.h"
+#include "MaterialLib/MPL/Property.h"
+
+namespace MaterialPropertyLib
+{
+class Medium;
+class Phase;
+class Component;
+
+/// Van Genuchten relative permeability function.
+///
+/// This property must be a medium property, it computes the permeability
+/// reduction due to saturation as function of capillary pressure.
+class RelPermVanGenuchten final : public Property
+{
+private:
+ double const _S_L_res;
+ double const _S_L_max;
+ double const _k_rel_min;
+ double const _m;
+ Medium* _medium = nullptr;
+
+public:
+ RelPermVanGenuchten(double const residual_liquid_saturation,
+ double const residual_gas_saturation,
+ double const min_relative_permeability_liquid,
+ double const exponent);
+ /// This method assigns a pointer to the material object that is the owner
+ /// of this property
+ void setScale(
+ std::variant<Medium*, Phase*, Component*> scale_pointer) override
+ {
+ if (std::holds_alternative<Medium*>(scale_pointer))
+ {
+ _medium = std::get<Medium*>(scale_pointer);
+ }
+ else
+ {
+ OGS_FATAL(
+ "The property 'RelativePermeabilityVanGenuchten' is "
+ "implemented on the 'media' scale only.");
+ }
+ }
+
+ /// Those methods override the base class implementations and
+ /// actually compute and set the property _values and _dValues.
+ PropertyDataType value(VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& pos,
+ double const t) const override;
+ PropertyDataType dValue(VariableArray const& variable_array,
+ Variable const variable,
+ ParameterLib::SpatialPosition const& pos,
+ double const t) const override;
+};
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/SaturationVanGenuchten.cpp b/MaterialLib/MPL/Properties/SaturationVanGenuchten.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7176b0bc65a226a4e3176408e5144eb4ae2b5227
--- /dev/null
+++ b/MaterialLib/MPL/Properties/SaturationVanGenuchten.cpp
@@ -0,0 +1,133 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "SaturationVanGenuchten.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include "MaterialLib/MPL/Medium.h"
+
+namespace MaterialPropertyLib
+{
+SaturationVanGenuchten::SaturationVanGenuchten(
+ double const residual_liquid_saturation,
+ double const residual_gas_saturation,
+ double const exponent,
+ double const entry_pressure)
+ : _S_L_res(residual_liquid_saturation),
+ _S_L_max(1. - residual_gas_saturation),
+ _m(exponent),
+ _p_b(entry_pressure)
+{
+ if (!(_m > 0 && _m < 1))
+ {
+ OGS_FATAL(
+ "The exponent value m = %g of van Genuchten saturation model, is "
+ "out of its range of (0, 1)",
+ _m);
+ }
+}
+
+PropertyDataType SaturationVanGenuchten::value(
+ VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& /*pos*/,
+ double const /*t*/) const
+{
+ const double p_cap = std::get<double>(
+ variable_array[static_cast<int>(Variable::capillary_pressure)]);
+
+ if (p_cap <= 0)
+ {
+ return _S_L_max;
+ }
+
+ double const p = p_cap / _p_b;
+ double const n = 1. / (1. - _m);
+ double const p_to_n = std::pow(p, n);
+
+ double const S_eff = std::pow(p_to_n + 1., -_m);
+ double const S = S_eff * _S_L_max - S_eff * _S_L_res + _S_L_res;
+ return std::clamp(S, _S_L_res, _S_L_max);
+}
+
+PropertyDataType SaturationVanGenuchten::dValue(
+ VariableArray const& variable_array, Variable const primary_variable,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/) const
+{
+ (void)primary_variable;
+ assert((primary_variable == Variable::capillary_pressure) &&
+ "SaturationVanGenuchten::dvalue is implemented for derivatives with "
+ "respect to capillary pressure only.");
+
+ const double p_cap = std::get<double>(
+ variable_array[static_cast<int>(Variable::capillary_pressure)]);
+
+ if (p_cap <= 0)
+ {
+ return 0;
+ }
+
+ double const p = p_cap / _p_b;
+ double const n = 1. / (1. - _m);
+ double const p_to_n = std::pow(p, n);
+
+ double const S_eff = std::pow(p_to_n + 1., -_m);
+ double const S = S_eff * _S_L_max - S_eff * _S_L_res + _S_L_res;
+
+ if (S < _S_L_res || S > _S_L_max)
+ {
+ return 0;
+ }
+
+ double const dS_eff_dp_cap = -_m * std::pow(p_cap / _p_b, n - 1) *
+ std::pow(1 + p_to_n, -1. - _m) /
+ (_p_b * (1. - _m));
+ return dS_eff_dp_cap * (_S_L_max - _S_L_res);
+}
+
+PropertyDataType SaturationVanGenuchten::d2Value(
+ VariableArray const& variable_array, Variable const primary_variable1,
+ Variable const primary_variable2,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/) const
+{
+ (void)primary_variable1;
+ (void)primary_variable2;
+ assert((primary_variable1 == Variable::capillary_pressure) &&
+ (primary_variable2 == Variable::capillary_pressure) &&
+ "SaturationVanGenuchten::d2Value is implemented for "
+ " derivatives with respect to capillary pressure only.");
+
+ const double p_cap = std::get<double>(
+ variable_array[static_cast<int>(Variable::capillary_pressure)]);
+
+ if (p_cap <= 0)
+ {
+ return 0;
+ }
+
+ double const p = p_cap / _p_b;
+ double const n = 1. / (1. - _m);
+ double const p_to_n = std::pow(p, n);
+
+ double const S_eff = std::pow(p_to_n + 1., -_m);
+ double const S = S_eff * _S_L_max - S_eff * _S_L_res + _S_L_res;
+
+ if (S < _S_L_res || S > _S_L_max)
+ {
+ return 0;
+ }
+
+ double const d2S_eff_dp_cap2 =
+ _m * p_to_n * std::pow(p_to_n + 1., -_m - 2.) * (p_to_n - _m) /
+ (p_cap * p_cap * (_m - 1.) * (_m - 1.));
+ return d2S_eff_dp_cap2 * (_S_L_max - _S_L_res);
+}
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/SaturationVanGenuchten.h b/MaterialLib/MPL/Properties/SaturationVanGenuchten.h
new file mode 100644
index 0000000000000000000000000000000000000000..2a792d3a7ce69cbe6f73b8404651dc9fc457eba2
--- /dev/null
+++ b/MaterialLib/MPL/Properties/SaturationVanGenuchten.h
@@ -0,0 +1,65 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ */
+#pragma once
+
+#include "MaterialLib/MPL/Property.h"
+
+namespace MaterialPropertyLib
+{
+class Medium;
+class Phase;
+class Component;
+
+/// The van Genuchten soil characteristics function.
+///
+/// This property must be a medium property, it computes the saturation of the
+/// wetting phase as function of capillary pressure.
+class SaturationVanGenuchten final : public Property
+{
+private:
+ Medium* _medium = nullptr;
+ double const _S_L_res;
+ double const _S_L_max;
+ double const _m;
+ double const _p_b;
+
+public:
+ SaturationVanGenuchten(double const residual_liquid_saturation,
+ double const residual_gas_saturation,
+ double const exponent,
+ double const entry_pressure);
+
+ void setScale(
+ std::variant<Medium*, Phase*, Component*> scale_pointer) override
+ {
+ if (!std::holds_alternative<Medium*>(scale_pointer))
+ {
+ OGS_FATAL(
+ "The property 'SaturationVanGenuchten' is implemented on the "
+ "'media' scale only.");
+ }
+ _medium = std::get<Medium*>(scale_pointer);
+ }
+
+ /// Those methods override the base class implementations and
+ /// actually compute and set the property _values and _dValues.
+ PropertyDataType value(VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& /*pos*/,
+ double const /*t*/) const override;
+ PropertyDataType dValue(VariableArray const& variable_array,
+ Variable const variable,
+ ParameterLib::SpatialPosition const& /*pos*/,
+ double const /*t*/) const override;
+ PropertyDataType d2Value(VariableArray const& variable_array,
+ Variable const variable1,
+ Variable const variable2,
+ ParameterLib::SpatialPosition const& /*pos*/,
+ double const /*t*/) const override;
+};
+} // namespace MaterialPropertyLib
diff --git a/Tests/MaterialLib/TestMPLRelPermVanGenuchten.cpp b/Tests/MaterialLib/TestMPLRelPermVanGenuchten.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cab51d5cc8f274464003ee85a4b6dccf2e0f1eff
--- /dev/null
+++ b/Tests/MaterialLib/TestMPLRelPermVanGenuchten.cpp
@@ -0,0 +1,72 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+#include <gtest/gtest.h>
+
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Properties/RelPermVanGenuchten.h"
+#include "TestMPL.h"
+#include "Tests/TestTools.h"
+
+namespace MPL = MaterialPropertyLib;
+
+TEST(MaterialPropertyLib, RelPermVanGenuchten)
+{
+ double const residual_liquid_saturation = 0.1;
+ double const residual_gas_saturation = 0.05;
+ double const exponent = 0.79;
+ double const min_relative_permeability_liquid = 1e-12;
+
+ MPL::Property const& permeability = MPL::RelPermVanGenuchten{
+ residual_liquid_saturation, residual_gas_saturation,
+ min_relative_permeability_liquid, exponent};
+
+ MPL::VariableArray variable_array;
+ ParameterLib::SpatialPosition const pos;
+ double const t = std::numeric_limits<double>::quiet_NaN();
+
+ double const S_0 = -0.1;
+ double const S_max = 1.1;
+ int const n_steps = 1000;
+ for (int i = 0; i <= n_steps; ++i)
+ {
+ double const S_L = S_0 + i * (S_max - S_0) / n_steps;
+ variable_array[static_cast<int>(MPL::Variable::liquid_saturation)] =
+ S_L;
+
+ double const k_rel =
+ permeability.template value<double>(variable_array, pos, t);
+ double const dk_rel = permeability.template dValue<double>(
+ variable_array, MPL::Variable::liquid_saturation, pos, t);
+
+ double const eps = 1e-8;
+ variable_array[static_cast<int>(MPL::Variable::liquid_saturation)] =
+ S_L - eps;
+ double const k_rel_minus =
+ permeability.template value<double>(variable_array, pos, t);
+ variable_array[static_cast<int>(MPL::Variable::liquid_saturation)] =
+ S_L + eps;
+ double const k_rel_plus =
+ permeability.template value<double>(variable_array, pos, t);
+
+ double const Dk_rel = (k_rel_plus - k_rel_minus) / 2 / eps;
+
+ if (S_L < 1 - residual_gas_saturation)
+ {
+ ASSERT_LE(std::abs(dk_rel - Dk_rel), 1e-7)
+ << "for saturation " << S_L << " and relative permeability "
+ << k_rel;
+ }
+ else
+ {
+ ASSERT_EQ(dk_rel, 0.) << "for saturation " << S_L
+ << " and relative permeability " << k_rel;
+ }
+ }
+}
diff --git a/Tests/MaterialLib/TestMPLSaturationVanGenuchten.cpp b/Tests/MaterialLib/TestMPLSaturationVanGenuchten.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f197e5f70083b6921dee10dabe6575d82884f251
--- /dev/null
+++ b/Tests/MaterialLib/TestMPLSaturationVanGenuchten.cpp
@@ -0,0 +1,69 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ */
+#include <gtest/gtest.h>
+
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Properties/SaturationVanGenuchten.h"
+#include "TestMPL.h"
+#include "Tests/TestTools.h"
+
+namespace MPL = MaterialPropertyLib;
+
+TEST(MaterialPropertyLib, SaturationVanGenuchten)
+{
+ double const residual_liquid_saturation = 0.1;
+ double const residual_gas_saturation = 0.05;
+ double const exponent = 0.79;
+ double const entry_pressure = 5000;
+
+ MPL::Property const& pressure_saturation = MPL::SaturationVanGenuchten{
+ residual_liquid_saturation, residual_gas_saturation, exponent,
+ entry_pressure};
+
+ MPL::VariableArray variable_array;
+ ParameterLib::SpatialPosition const pos;
+ double const t = std::numeric_limits<double>::quiet_NaN();
+
+ double const p_0 = -1e6;
+ double const p_max = 10000;
+ int const n_steps = 10000;
+ for (int i = 0; i <= n_steps; ++i)
+ {
+ double const p_L = p_0 + i * (p_max - p_0) / n_steps;
+ variable_array[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ -p_L;
+
+ double const S =
+ pressure_saturation.template value<double>(variable_array, pos, t);
+ double const dS = pressure_saturation.template dValue<double>(
+ variable_array, MPL::Variable::capillary_pressure, pos, t);
+ double const dS2 = pressure_saturation.template d2Value<double>(
+ variable_array, MPL::Variable::capillary_pressure,
+ MPL::Variable::capillary_pressure, pos, t);
+
+ double const eps = 1e-1;
+ variable_array[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ -p_L - eps;
+ double const S_minus =
+ pressure_saturation.template value<double>(variable_array, pos, t);
+ variable_array[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ -p_L + eps;
+ double const S_plus =
+ pressure_saturation.template value<double>(variable_array, pos, t);
+
+ double const DS = (S_plus - S_minus) / 2 / eps;
+ double const DS2 = (S_plus - 2 * S + S_minus) / (eps * eps);
+
+ ASSERT_LE(std::abs(dS - DS), 1e-9)
+ << "for capillary pressure " << -p_L << " and saturation " << S;
+ ASSERT_LE(std::abs(dS2 - DS2), 1e-9)
+ << "for capillary pressure " << -p_L << " and saturation " << S;
+ }
+}