diff --git a/Documentation/ProjectFile/properties/property/SaturationExponential/c_SaturationExponential.md b/Documentation/ProjectFile/properties/property/SaturationExponential/c_SaturationExponential.md
new file mode 100644
index 0000000000000000000000000000000000000000..83ef80845282aff21a3c637df424d9983b54f811
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/SaturationExponential/c_SaturationExponential.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::SaturationExponential
diff --git a/Documentation/ProjectFile/properties/property/SaturationExponential/t_exponent.md b/Documentation/ProjectFile/properties/property/SaturationExponential/t_exponent.md
new file mode 100644
index 0000000000000000000000000000000000000000..baa7d302274f8a601376532e568ca9168413820f
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/SaturationExponential/t_exponent.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::SaturationExponential::exponent_
diff --git a/Documentation/ProjectFile/properties/property/SaturationExponential/t_maximum_capillary_pressure.md b/Documentation/ProjectFile/properties/property/SaturationExponential/t_maximum_capillary_pressure.md
new file mode 100644
index 0000000000000000000000000000000000000000..52d2d511621d6a12a6d7066e129a396fed49a31d
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/SaturationExponential/t_maximum_capillary_pressure.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::SaturationExponential::p_cap_max_
diff --git a/Documentation/ProjectFile/properties/property/SaturationExponential/t_residual_gas_saturation.md b/Documentation/ProjectFile/properties/property/SaturationExponential/t_residual_gas_saturation.md
new file mode 100644
index 0000000000000000000000000000000000000000..0b56bfad36cc6da26f73a9dcf10a7e2f58226dd8
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/SaturationExponential/t_residual_gas_saturation.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::SaturationExponential::residual_gas_saturation_
diff --git a/Documentation/ProjectFile/properties/property/SaturationExponential/t_residual_liquid_saturation.md b/Documentation/ProjectFile/properties/property/SaturationExponential/t_residual_liquid_saturation.md
new file mode 100644
index 0000000000000000000000000000000000000000..ce79c90316d14b49e7c2790e812ab67db34d0ffb
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/SaturationExponential/t_residual_liquid_saturation.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::SaturationExponential::residual_liquid_saturation_
diff --git a/MaterialLib/MPL/CreateProperty.cpp b/MaterialLib/MPL/CreateProperty.cpp
index 789572824f6553bc25a17c17476b3b4592ba9500..76975329ab4c2a1a5cb11e02cd207aa46bd68343 100644
--- a/MaterialLib/MPL/CreateProperty.cpp
+++ b/MaterialLib/MPL/CreateProperty.cpp
@@ -157,6 +157,10 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
return createRelPermLiakopoulos(config);
}
+ if (boost::iequals(property_type, "SaturationExponential"))
+ {
+ return createSaturationExponential(config);
+ }
if (boost::iequals(property_type, "SaturationVanGenuchten"))
{
return createSaturationVanGenuchten(config);
diff --git a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationExponential.cpp b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationExponential.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8909044823df4f543a0b08d1dff6e666dbc6ad97
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationExponential.cpp
@@ -0,0 +1,45 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 "SaturationExponential.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<SaturationExponential> createSaturationExponential(
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{properties__property__type}
+ config.checkConfigParameter("type", "SaturationExponential");
+
+ // Second access for storage.
+ //! \ogs_file_param{properties__property__name}
+ auto property_name = config.peekConfigParameter<std::string>("name");
+
+ DBUG("Create SaturationExponential medium property {:s}.", property_name);
+
+ auto const residual_liquid_saturation =
+ //! \ogs_file_param{properties__property__SaturationExponential__residual_liquid_saturation}
+ config.getConfigParameter<double>("residual_liquid_saturation");
+ auto const residual_gas_saturation =
+ //! \ogs_file_param{properties__property__SaturationExponential__residual_gas_saturation}
+ config.getConfigParameter<double>("residual_gas_saturation");
+ auto const p_cap_max =
+ //! \ogs_file_param{properties__property__SaturationExponential__maximum_capillary_pressure}
+ config.getConfigParameter<double>("maximum_capillary_pressure");
+ auto const exponent =
+ //! \ogs_file_param{properties__property__SaturationExponential__exponent}
+ config.getConfigParameter<double>("exponent");
+
+ return std::make_unique<SaturationExponential>(
+ std::move(property_name), residual_liquid_saturation,
+ residual_gas_saturation, p_cap_max, exponent);
+}
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationExponential.h b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationExponential.h
new file mode 100644
index 0000000000000000000000000000000000000000..2e00abdee2f6b594d97afb4ec0b77c702d08d811
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/CreateSaturationExponential.h
@@ -0,0 +1,29 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 SaturationExponential;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<SaturationExponential> createSaturationExponential(
+ BaseLib::ConfigTree const& config);
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.cpp b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f9bbe5a9f7fd6aee827d8b5ff0898e2f337d51d7
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.cpp
@@ -0,0 +1,94 @@
+/**
+ * \file
+ * \author Norbert Grunwald
+ * \date 27.06.2018
+ * \brief
+ *
+ * \copyright
+ * Copyright (c) 2012-2021, 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 "SaturationExponential.h"
+
+#include <algorithm>
+#include <cmath>
+
+#include "MaterialLib/MPL/Medium.h"
+#include "MathLib/MathTools.h"
+
+namespace MaterialPropertyLib
+{
+SaturationExponential::SaturationExponential(
+ std::string name,
+ const double residual_liquid_saturation,
+ const double residual_gas_saturation,
+ const double p_cap_max,
+ const double exponent)
+ : residual_liquid_saturation_(residual_liquid_saturation),
+ residual_gas_saturation_(residual_gas_saturation),
+ p_cap_max_(p_cap_max),
+ exponent_(exponent)
+{
+ name_ = std::move(name);
+ if (p_cap_max_ <= 0.)
+ {
+ OGS_FATAL(
+ "Reference capillary pressure must be positive in "
+ "MPL::SaturationExponential.");
+ }
+};
+
+PropertyDataType SaturationExponential::value(
+ VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+ double const /*dt*/) const
+{
+ const double p_cap = std::get<double>(
+ variable_array[static_cast<int>(Variable::capillary_pressure)]);
+
+ const double S_res = residual_liquid_saturation_;
+ const double S_max = 1. - residual_gas_saturation_;
+
+ const double pc = std::clamp(p_cap, 0., p_cap_max_);
+ const double S_e = 1. - std::pow(pc / p_cap_max_, exponent_);
+ return S_e * (S_max - S_res) + S_res;
+}
+
+PropertyDataType SaturationExponential::dValue(
+ VariableArray const& variable_array, Variable const primary_variable,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+ double const /*dt*/) const
+{
+ (void)primary_variable;
+ assert((primary_variable == Variable::capillary_pressure) &&
+ "SaturationExponential::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)]);
+
+ const double S_res = residual_liquid_saturation_;
+ const double S_max = 1. - residual_gas_saturation_;
+
+ if ((p_cap > p_cap_max_) || (p_cap <= 0.))
+ {
+ return 0.;
+ }
+ return (exponent_ / p_cap) * (S_res - S_max) *
+ std::pow(p_cap / p_cap_max_, exponent_);
+}
+
+PropertyDataType SaturationExponential::d2Value(
+ VariableArray const& /*variable_array*/,
+ Variable const /*primary_variable1*/, Variable const /*primary_variable2*/,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+ double const /*dt*/) const
+{
+ OGS_FATAL("SaturationExponential::d2Value() is not implemented.");
+}
+
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.h b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.h
new file mode 100644
index 0000000000000000000000000000000000000000..52b02c6f46205f6b511b3f3a6e8766c4328182fb
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.h
@@ -0,0 +1,82 @@
+/**
+ * \file
+ * \author Norbert Grunwald
+ * \date 27.06.2018
+ * \brief
+ *
+ * \copyright
+ * Copyright (c) 2012-2021, 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
+{
+/**
+ * \brief A simplistic soil characteristics function.
+ * \details This property must be a medium property. It
+ * computes the saturation of the wetting phase as function
+ * of capillary pressure.
+ *
+ * The capillary pressure-saturation relationship given by
+ * \f$S_{eff}=1-\left(\frac{p_{c}}{p_{c}^{max}}\right)^{\lambda}\f$,
+ * where
+ * \f$\lambda\f$ is an exponent,
+ * \f$p_{c}\f$ is capillary pressure,
+ * \f$p_{c}^{max}\f$ is the maximum capillary pressure at which \f$S_{eff}=0\f$.
+ *
+ * This property can mainly be used for testing. If the exponent is set to 1,
+ * the characteristic curve shows a linear dependence.
+ *
+ *
+ * Reference capillary pressure at which \f$S_{eff}=0\f$.
+ */
+class SaturationExponential final : public Property
+{
+private:
+ const double
+ residual_liquid_saturation_; ///< Residual saturation of the gas phase.
+ const double
+ residual_gas_saturation_; ///< Residual saturation of the liquid phase.
+ const double p_cap_max_; ///< Maximum capillary pressure at which
+ /// effective saturation is zero.
+ const double exponent_; ///< Exponent to govern the shape of the curve.
+
+public:
+ SaturationExponential(std::string name,
+ const double residual_liquid_saturation,
+ const double residual_gas_saturation,
+ const double p_cap_max, const double exponent);
+
+ void checkScale() const override
+ {
+ if (!std::holds_alternative<Medium*>(scale_))
+ {
+ OGS_FATAL(
+ "The property 'SaturationExponential' is implemented on the "
+ "'media' scale only.");
+ }
+ }
+
+ PropertyDataType value(VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& /*pos*/,
+ double const /*t*/,
+ double const /*dt*/) const override;
+ PropertyDataType dValue(VariableArray const& variable_array,
+ Variable const primary_variable,
+ ParameterLib::SpatialPosition const& /*pos*/,
+ double const /*t*/,
+ double const /*dt*/) const override;
+ PropertyDataType d2Value(VariableArray const& variable_array,
+ Variable const /*primary_variable1*/,
+ Variable const /*primary_variable2*/,
+ ParameterLib::SpatialPosition const& /*pos*/,
+ double const /*t*/,
+ double const /*dt*/) const override;
+};
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CreateProperties.h b/MaterialLib/MPL/Properties/CreateProperties.h
index 705477492ac97a86837bb3a80d3eb19b505d8311..5cf7b9628bcbff60ba870dff96ff00057f3f6774 100644
--- a/MaterialLib/MPL/Properties/CreateProperties.h
+++ b/MaterialLib/MPL/Properties/CreateProperties.h
@@ -15,6 +15,7 @@
#include "CapillaryPressureSaturation/CreateCapillaryPressureVanGenuchten.h"
#include "CapillaryPressureSaturation/CreateSaturationBrooksCorey.h"
#include "CapillaryPressureSaturation/CreateSaturationLiakopoulos.h"
+#include "CapillaryPressureSaturation/CreateSaturationExponential.h"
#include "CapillaryPressureSaturation/CreateSaturationVanGenuchten.h"
#include "CreateAverageMolarMass.h"
#include "CreateBishopsPowerLaw.h"
diff --git a/MaterialLib/MPL/Properties/Properties.h b/MaterialLib/MPL/Properties/Properties.h
index 466e88f407384aeb28d64b1c62e8277949f3b74f..4af58ab9e4714442e3320a61990844942de7d53d 100644
--- a/MaterialLib/MPL/Properties/Properties.h
+++ b/MaterialLib/MPL/Properties/Properties.h
@@ -16,6 +16,7 @@
#include "BishopsSaturationCutoff.h"
#include "CapillaryPressureSaturation/SaturationBrooksCorey.h"
#include "CapillaryPressureSaturation/SaturationLiakopoulos.h"
+#include "CapillaryPressureSaturation/SaturationExponential.h"
#include "CapillaryPressureSaturation/SaturationVanGenuchten.h"
#include "ClausiusClapeyron.h"
#include "Constant.h"
diff --git a/Tests/MaterialLib/TestMPLSaturationExponential.cpp b/Tests/MaterialLib/TestMPLSaturationExponential.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e1af368eda17f0d4fd616917f398223a6c4dd4c5
--- /dev/null
+++ b/Tests/MaterialLib/TestMPLSaturationExponential.cpp
@@ -0,0 +1,82 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2021, 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 <iomanip>
+
+#include "MaterialLib/MPL/Medium.h"
+#include "MaterialLib/MPL/Properties/CapillaryPressureSaturation/SaturationExponential.h"
+#include "TestMPL.h"
+#include "Tests/TestTools.h"
+
+namespace MPL = MaterialPropertyLib;
+
+TEST(MaterialPropertyLib, SaturationExponential)
+{
+ double const residual_liquid_saturation = 0.1;
+ double const residual_gas_saturation = 0.05;
+ double const exponent = 1.5;
+ double const p_cap_max = 30000;
+
+ MPL::Property const& saturation = MPL::SaturationExponential{
+ "saturation", residual_liquid_saturation, residual_gas_saturation,
+ p_cap_max, exponent};
+
+ MPL::VariableArray variable_array;
+ ParameterLib::SpatialPosition const pos;
+ double const t = std::numeric_limits<double>::quiet_NaN();
+ double const dt = std::numeric_limits<double>::quiet_NaN();
+
+ double const p_0 = -p_cap_max;
+ double const p_max = 1.5 * p_cap_max;
+ int const n_steps = 9999;
+ for (int i = 0; i <= n_steps; ++i)
+ {
+ double const p_cap = p_0 + i * (p_max - p_0) / n_steps;
+ variable_array[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ p_cap;
+
+ double const s_res = residual_liquid_saturation;
+ double const s_max = 1. - residual_gas_saturation;
+ double const s_eff =
+ 1. -
+ std::pow(std::clamp(p_cap, 0., p_cap_max) / p_cap_max, exponent);
+
+ double s_ref = s_eff * (s_max - s_res) + s_res;
+
+ double const S =
+ saturation.template value<double>(variable_array, pos, t, dt);
+ double const dS = saturation.template dValue<double>(
+ variable_array, MPL::Variable::capillary_pressure, pos, t, dt);
+
+ double const eps = 1e-2;
+ variable_array[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ p_cap - eps;
+
+ double const S_minus =
+ saturation.template value<double>(variable_array, pos, t, dt);
+
+ variable_array[static_cast<int>(MPL::Variable::capillary_pressure)] =
+ p_cap + eps;
+
+ double const S_plus =
+ saturation.template value<double>(variable_array, pos, t, dt);
+
+ double const dS_ref = p_cap > 0 ? (S_plus - S_minus) / 2 / eps : 0.;
+
+ ASSERT_LE(std::abs(s_ref - S), 1e-9)
+ << std::setprecision(16) << "with: \ncapillary pressure: " << p_cap
+ << "\nsaturation: " << S << "\nsaturation_ref = " << s_ref
+ << "\nat point: " << i;
+ ASSERT_LE(std::abs(dS - dS_ref), 1e-9)
+ << std::setprecision(16) << "with: \np_cap: " << p_cap
+ << "\nS: " << S << "\nS_ref = " << s_ref << "\ndS = " << dS
+ << "\ndS_ref = " << dS_ref << "\nat point: " << i;
+ }
+}