diff --git a/Documentation/ProjectFile/properties/property/t_WaterThermalConductivityIAPWS.md b/Documentation/ProjectFile/properties/property/t_WaterThermalConductivityIAPWS.md
new file mode 100644
index 0000000000000000000000000000000000000000..bbf198ef5cb02631497460d72eb4ec1d238da903
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/t_WaterThermalConductivityIAPWS.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::WaterThermalConductivityIAPWS
diff --git a/MaterialLib/MPL/CreateProperty.cpp b/MaterialLib/MPL/CreateProperty.cpp
index c268ee7ad375902fcd9d63eb9a20017ade4d287d..5d723f4e253e85d7f8008af1ec6401f13289a75e 100644
--- a/MaterialLib/MPL/CreateProperty.cpp
+++ b/MaterialLib/MPL/CreateProperty.cpp
@@ -300,16 +300,21 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
{
return createVolumeFractionAverage(config);
}
+
if (boost::iequals(property_type, "WaterViscosityIAPWS"))
{
return createWaterViscosityIAPWS(config);
}
+
if (boost::iequals(property_type, "LiquidViscosityVogels"))
{
return createLiquidViscosityVogels(config);
}
-
+ if (boost::iequals(property_type, "WaterThermalConductivityIAPWS"))
+ {
+ return createWaterThermalConductivityIAPWS(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",
diff --git a/MaterialLib/MPL/Properties/CreateProperties.h b/MaterialLib/MPL/Properties/CreateProperties.h
index c3522b4fb1aefeea6762482ac73623ff73397f6c..8f5a19a5a11ab81f8fac56cfc7ca8ea2aba06854 100644
--- a/MaterialLib/MPL/Properties/CreateProperties.h
+++ b/MaterialLib/MPL/Properties/CreateProperties.h
@@ -65,3 +65,4 @@
#include "VapourDiffusion/CreateVapourDiffusionPMQ.h"
#include "Viscosity/CreateWaterViscosityIAPWS.h"
#include "Viscosity/CreateLiquidViscosityVogels.h"
+#include "ThermalConductivity/CreateWaterThermalConductivityIAPWS.h"
diff --git a/MaterialLib/MPL/Properties/ThermalConductivity/CreateWaterThermalConductivityIAPWS.cpp b/MaterialLib/MPL/Properties/ThermalConductivity/CreateWaterThermalConductivityIAPWS.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..467f22153feca2b55abf2b439aa83fbcf3e41874
--- /dev/null
+++ b/MaterialLib/MPL/Properties/ThermalConductivity/CreateWaterThermalConductivityIAPWS.cpp
@@ -0,0 +1,34 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * Created on March 4, 2021, 4:38 PM
+ */
+
+#include "CreateWaterThermalConductivityIAPWS.h"
+
+#include "BaseLib/ConfigTree.h"
+#include "MaterialLib/MPL/Property.h"
+#include "WaterThermalConductivityIAPWS.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<Property> createWaterThermalConductivityIAPWS(
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{properties__property__type}
+ config.checkConfigParameter("type", "WaterThermalConductivityIAPWS");
+ DBUG("Create WaterThermalConductivityIAPWS phase property");
+
+ // Second access for storage.
+ //! \ogs_file_param{properties__property__name}
+ auto property_name = config.peekConfigParameter<std::string>("name");
+
+ //! \ogs_file_param_special{properties__property__WaterThermalConductivityIAPWS}
+ return std::make_unique<WaterThermalConductivityIAPWS>(std::move(property_name));
+}
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/ThermalConductivity/CreateWaterThermalConductivityIAPWS.h b/MaterialLib/MPL/Properties/ThermalConductivity/CreateWaterThermalConductivityIAPWS.h
new file mode 100644
index 0000000000000000000000000000000000000000..f5bbbdfa2d7f5d36313fada78626119bc6a3adf1
--- /dev/null
+++ b/MaterialLib/MPL/Properties/ThermalConductivity/CreateWaterThermalConductivityIAPWS.h
@@ -0,0 +1,26 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * Created on March 4, 2021, 4:38 PM
+ */
+
+#pragma once
+
+#include <memory>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+class Property;
+std::unique_ptr<Property> createWaterThermalConductivityIAPWS(
+ BaseLib::ConfigTree const& config);
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/ThermalConductivity/WaterThermalConductivityIAPWS.cpp b/MaterialLib/MPL/Properties/ThermalConductivity/WaterThermalConductivityIAPWS.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9dab62ce530ecda494bf48cd51afb0bf526b24bd
--- /dev/null
+++ b/MaterialLib/MPL/Properties/ThermalConductivity/WaterThermalConductivityIAPWS.cpp
@@ -0,0 +1,201 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * Created on March 4, 2021, 3:05 PM
+ */
+
+#include "WaterThermalConductivityIAPWS.h"
+
+#include <cmath>
+
+#include "BaseLib/Error.h"
+#include "MaterialLib/MPL/Medium.h"
+
+namespace MaterialPropertyLib
+{
+ // Li and Lij are coefficients from Tables 1 and 2 (Daucik and Dooley, 2011)
+static const double Li[5] = {2.443221e-3, 1.323095e-2, 6.770357e-3, -3.454586e-3, 4.096266e-4};
+static const double Lij[5][6] = {
+ {1.60397357, -0.646013523, 0.111443906, 0.102997357, -0.0504123634, 0.00609859258},
+ {2.33771842, -2.78843778, 1.53616167, -0.463045512, 0.0832827019, -0.00719201245},
+ {2.19650529, -4.54580785, 3.55777244, -1.40944978, 0.275418278, -0.0205938816},
+ {-1.21051378, 1.60812989, -0.621178141, 0.0716373224, 0, 0},
+ {-2.7203370, 4.57586331, -3.18369245, 1.1168348, -0.19268305, 0.012913842}};
+
+static double computeBarLambda0Factor(const double barT);
+
+static std::array<double, 5> computeSeriesFactorTForLambda1(const double barT);
+static std::array<double, 6> computeSeriesFactorRhoForLambda1(const double bar_rho);
+static double computeBarLambda1Factor(
+ const std::array<double, 5>& series_factorT,
+ const std::array<double, 6>& series_factorRho);
+
+static double computedBarLambda_dbarT(const double barT, double bar_rho);
+static double computedBarLambda_dbarRho(const double barT, double bar_rho);
+
+
+PropertyDataType WaterThermalConductivityIAPWS::value(
+ VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+ double const /*dt*/) const
+{
+ const double bar_T = variable_array.temperature / ref_T_;
+ const double bar_rho = variable_array.density / ref_rho_;
+
+ const double lambda0 = std::sqrt(bar_T) / computeBarLambda0Factor(bar_T);
+
+ const auto& series_factorT = computeSeriesFactorTForLambda1(bar_T);
+ const auto& series_factorRho = computeSeriesFactorRhoForLambda1(bar_rho);
+ const double lambda1 = std::exp(
+ bar_rho * computeBarLambda1Factor(series_factorT, series_factorRho));
+
+ return lambda0 * lambda1 * ref_lambda_;
+}
+
+PropertyDataType WaterThermalConductivityIAPWS::dValue(
+ VariableArray const& variable_array, Variable const variable,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+ double const /*dt*/) const
+{
+ const double bar_T = variable_array.temperature / ref_T_;
+ const double bar_rho = variable_array.density / ref_rho_;
+
+ switch (variable)
+ {
+ case Variable::temperature:
+ return ref_lambda_ * computedBarLambda_dbarT(bar_T, bar_rho) / ref_T_;
+ case Variable::density:
+ return ref_lambda_ * computedBarLambda_dbarRho(bar_T, bar_rho) / ref_rho_;
+ default:
+ OGS_FATAL(
+ "WaterThermalConductivityIAPWS::dValue is implemented for "
+ "derivatives with respect to temperature and liquid density only.");
+ }
+}
+
+double computeBarLambda0Factor(const double barT)
+{
+ double sum_val = 0.;
+ double barT_i = 1.;
+ for (double value : Li)
+ {
+ sum_val += (value / barT_i);
+ barT_i *= barT;
+ }
+ return sum_val;
+}
+
+std::array<double, 5> computeSeriesFactorTForLambda1(const double barT)
+{
+ std::array<double, 5> series_factorT;
+ series_factorT[0] = 1.;
+ const double barT_fac = 1 / barT - 1.0;
+ for (int i = 1; i < 5; i++)
+ {
+ series_factorT[i] = series_factorT[i - 1] * barT_fac;
+ }
+
+ return series_factorT;
+}
+
+std::array<double, 6> computeSeriesFactorRhoForLambda1(const double bar_rho)
+{
+ std::array<double, 6> series_factorRho;
+ series_factorRho[0] = 1.;
+ for (int i = 1; i < 6; i++)
+ {
+ series_factorRho[i] = series_factorRho[i - 1] * (bar_rho - 1.0);
+ }
+
+ return series_factorRho;
+}
+
+double computeBarLambda1Factor(const std::array<double, 5>& series_factorT,
+ const std::array<double, 6>& series_factorRho)
+{
+ double sum_val = 0.;
+ for (int i = 0; i < 5; i++)
+ {
+ double sum_val_j = 0;
+ for (int j = 0; j < 6; j++)
+ {
+ sum_val_j += Lij[i][j] * series_factorRho[j];
+ }
+ sum_val += series_factorT[i] * sum_val_j;
+ }
+
+ return sum_val;
+}
+
+double computedBarLambda_dbarT(const double barT, double bar_rho)
+{
+ const double lambda0_factor = computeBarLambda0Factor(barT);
+ const double sqrt_barT = std::sqrt(barT);
+
+ double dlambda0_factor_dbarT = 0.0;
+ double barT_i = barT * barT;
+ for (int i = 1; i < 5; i++)
+ {
+ dlambda0_factor_dbarT -= static_cast<double>(i) * (Li[i] / barT_i);
+ barT_i *= barT;
+ }
+
+ const double dbar_lambda0_dbarT =
+ 0.5 / (lambda0_factor * sqrt_barT) -
+ sqrt_barT * dlambda0_factor_dbarT / (lambda0_factor * lambda0_factor);
+
+ const auto& series_factorT = computeSeriesFactorTForLambda1(barT);
+ const auto& series_factorRho = computeSeriesFactorRhoForLambda1(bar_rho);
+
+ double dlambda1_factor_dbarT = 0.0;
+ for (int i = 1; i < 5; i++)
+ {
+ double sum_val_j = 0;
+ for (int j = 0; j < 6; j++)
+ {
+ sum_val_j += Lij[i][j] * series_factorRho[j];
+ }
+ dlambda1_factor_dbarT -= static_cast<double>(i) * series_factorT[i - 1] *
+ sum_val_j / (barT * barT);
+ }
+
+ const double lambda1_factor =
+ computeBarLambda1Factor(series_factorT, series_factorRho);
+ const double dbar_lambda1_dbarT =
+ bar_rho * std::exp(bar_rho * lambda1_factor) * dlambda1_factor_dbarT;
+
+ return dbar_lambda0_dbarT * std::exp(bar_rho * lambda1_factor) +
+ dbar_lambda1_dbarT * sqrt_barT / lambda0_factor;
+}
+
+double computedBarLambda_dbarRho(const double barT, double bar_rho)
+{
+ const auto& series_factorT = computeSeriesFactorTForLambda1(barT);
+ const auto& series_factorRho = computeSeriesFactorRhoForLambda1(bar_rho);
+
+ double dlambda1_factor_dbar_rho = 0.0;
+ for (int i = 0; i < 5; i++)
+ {
+ double sum_val_j = 0;
+ for (int j = 1; j < 6; j++)
+ {
+ sum_val_j +=
+ static_cast<double>(j) * Lij[i][j] * series_factorRho[j - 1];
+ }
+ dlambda1_factor_dbar_rho += series_factorT[i] * sum_val_j;
+ }
+
+ const double lambda0 = std::sqrt(barT) / computeBarLambda0Factor(barT);
+
+ const double lambda1_factor =
+ computeBarLambda1Factor(series_factorT, series_factorRho);
+ return lambda0 * std::exp(bar_rho * lambda1_factor) *
+ (lambda1_factor + bar_rho * dlambda1_factor_dbar_rho);
+}
+
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/ThermalConductivity/WaterThermalConductivityIAPWS.h b/MaterialLib/MPL/Properties/ThermalConductivity/WaterThermalConductivityIAPWS.h
new file mode 100644
index 0000000000000000000000000000000000000000..1324fcb87896238eb0807903b52174f02eb6940f
--- /dev/null
+++ b/MaterialLib/MPL/Properties/ThermalConductivity/WaterThermalConductivityIAPWS.h
@@ -0,0 +1,63 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
+ * Distributed under a Modified BSD License.
+ * See accompanying file LICENSE.txt or
+ * http://www.opengeosys.org/project/license
+ *
+ * Created on March 4, 2021, 3:05 PM
+ */
+
+#pragma once
+
+#include "MaterialLib/MPL/Property.h"
+
+namespace MaterialPropertyLib
+{
+class Phase;
+
+/**
+ * \brief A class for thermal conductivity model that is defined by
+ * The International Association for the Properties of Water and Steam
+ * <a href="http://www.iapws.org/relguide/ThCond.pdf">IAPWS</a>
+ * (File accessed at 13.01.2023) - (Daucik and Dooley, 2011)
+ *
+ * With the definition, the thermal conductivity is a function of temperature and
+ * water density
+ *
+ * \attention The critical enhancement, \f$\bar{\lambda}_2\f$, is not considered.
+ * For information on region of significance and the significance,
+ * please see Figure 2 from the document linked in the upper paragraph.
+ */
+class WaterThermalConductivityIAPWS final : public Property
+{
+public:
+ explicit WaterThermalConductivityIAPWS(std::string name) { name_ = std::move(name); }
+ void checkScale() const override
+ {
+ if (!std::holds_alternative<Phase*>(scale_))
+ {
+ OGS_FATAL(
+ "The property 'WaterThermalConductivityIAPWS' is "
+ "implemented on the 'Phase' scale only.");
+ }
+ }
+
+ /// \return The water thermal conductivity.
+ PropertyDataType value(VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& pos,
+ double const t, double const dt) const override;
+ /// \return The derivative of water thermal conductivity with respect to
+ /// temperature or phase (water) pressure.
+ PropertyDataType dValue(VariableArray const& variable_array,
+ Variable const variable,
+ ParameterLib::SpatialPosition const& pos,
+ double const t, double const dt) const override;
+private:
+ static constexpr double ref_T_ = 647.096; ///< reference temperature in K
+ static constexpr double ref_rho_ = 322.0; ///< reference density in `kg/m^3`
+ static constexpr double ref_lambda_ = 1.0e-3; ///< reference thermal conductivity in `W.K^-1.m^-1`
+
+};
+} // namespace MaterialPropertyLib
diff --git a/Tests/MaterialLib/TestMPLWaterThermalConductivity.cpp b/Tests/MaterialLib/TestMPLWaterThermalConductivity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..add5fcaf9db44a28ceb933c15123ad74258b30a8
--- /dev/null
+++ b/Tests/MaterialLib/TestMPLWaterThermalConductivity.cpp
@@ -0,0 +1,78 @@
+/**
+ * \file
+ * \brief Test thermal conductivity of water IAPWS97 model
+ *
+ * \copyright
+ * Copyright (c) 2012-2023, 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 <memory>
+
+#include "MaterialLib/MPL/Properties/ThermalConductivity/CreateWaterThermalConductivityIAPWS.h"
+#include "TestMPL.h"
+
+TEST(Material, checkWaterThermalConductivityIAPWS_)
+{
+ const char xml[] =
+ "<property>"
+ " <name>thermal_conductivity</name>"
+ " <type>WaterThermalConductivityIAPWS</type>"
+ "</property>";
+
+ std::unique_ptr<MaterialPropertyLib::Property> const property_ptr =
+ Tests::createTestProperty(
+ xml, MaterialPropertyLib::createWaterThermalConductivityIAPWS);
+ MaterialPropertyLib::Property const& property = *property_ptr;
+
+ MaterialPropertyLib::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();
+
+ // Test data provided on http://www.iapws.org/relguide/ThCond.pdf Table 4
+ const double T[] = {298.15, 298.15, 298.15, 873.15};
+ const double rho[] = {0, 998, 1200, 0};
+
+ const double expected_lambda[] = {18.4341883e-3, 607.712868e-3, 799.038144e-3, 79.1034659e-3};
+
+ const double perturbation = 1.e-8;
+ for (int i = 0; i < 4; i++)
+ {
+ // Test lambda
+ variable_array.temperature = T[i];
+ variable_array.density = rho[i];
+ const double lambda = property.template value<double>(
+ variable_array, pos, t, dt);
+ ASSERT_NEAR(expected_lambda[i], lambda, 1.e-9);
+
+ const double dlambda_dT = property.template dValue<double>(
+ variable_array, MaterialPropertyLib::Variable::temperature, pos,
+ t, dt);
+ const double dlambda_drho = property.template dValue<double>(
+ variable_array, MaterialPropertyLib::Variable::density, pos,
+ t, dt);
+
+ // Test dlambda/dT
+ variable_array.temperature =
+ T[i] + perturbation;
+ double lambda1 = property.template value<double>(
+ variable_array, pos, t, dt);
+ ASSERT_NEAR((lambda1 - lambda) / perturbation, dlambda_dT, 8.e-6);
+
+ // Test dlambda/drho
+ variable_array.temperature = T[i];
+ variable_array.density =
+ rho[i] + perturbation;
+ lambda1 = property.template value<double>(
+ variable_array, pos, t, dt);
+
+ ASSERT_NEAR((lambda1 - lambda) / perturbation, dlambda_drho, 8.e-6);
+ }
+}