diff --git a/Tests/MaterialLib/TestMPLWaterVapourDensity.cpp b/Tests/MaterialLib/TestMPLWaterVapourDensity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..113c1685d1f0b47be567727931a74043f64aef19
--- /dev/null
+++ b/Tests/MaterialLib/TestMPLWaterVapourDensity.cpp
@@ -0,0 +1,149 @@
+/**
+ * \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
+ *
+ * Created on March 4, 2021, 5:23 PM
+ */
+#include <gtest/gtest.h>
+
+#include <array>
+#include <cmath>
+#include <limits>
+
+#include "MaterialLib/MPL/Properties/Density/CreateWaterVapourDensity.h"
+#include "TestMPL.h"
+
+TEST(MaterialPropertyLib, WaterVapourDensity)
+{
+ char const xml[] =
+ "<property>"
+ " <name>vapour_density</name>"
+ " <type>WaterVapourDensity</type>"
+ "</property>";
+
+ std::unique_ptr<MaterialPropertyLib::Property> const property_ptr =
+ Tests::createTestProperty(
+ xml, MaterialPropertyLib::createWaterVapourDensity);
+ MaterialPropertyLib::Property const& property = *property_ptr;
+
+ double const T = 390.0;
+ double const p = -1.e+6;
+ double const rho_w = 1000.0;
+
+ 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();
+ variable_array[static_cast<int>(
+ MaterialPropertyLib::Variable::temperature)] = T;
+ variable_array[static_cast<int>(
+ MaterialPropertyLib::Variable::phase_pressure)] = p;
+ variable_array[static_cast<int>(MaterialPropertyLib::Variable::density)] =
+ rho_w;
+
+ // The derivative of the water vapour with respect of temperature
+ {
+ std::array const temperatures = {273.0, 293.0, 393.0, 420.0, 500.0};
+ std::array const rho_vw_expected = {4.875989e-03, 1.692871e-02,
+ 1.276865, 2.882635, 1.920386e+01};
+
+ for (std::size_t i = 0; i < temperatures.size(); ++i)
+ {
+ double const T_i = temperatures[i];
+ variable_array[static_cast<int>(
+ MaterialPropertyLib::Variable::temperature)] = T_i;
+
+ double const rho_vw =
+ property.template value<double>(variable_array, pos, t, dt);
+
+ ASSERT_LE(std::fabs(rho_vw_expected[i] - rho_vw), 5e-6)
+ << "for expected water vapour density " << rho_vw_expected[i]
+ << " and for computed water vapour density " << rho_vw;
+
+ double const dT = 1.0e-4;
+ variable_array[static_cast<int>(
+ MaterialPropertyLib::Variable::temperature)] = T_i - dT;
+ double const rho_vw0 =
+ property.template value<double>(variable_array, pos, t, dt);
+
+ variable_array[static_cast<int>(
+ MaterialPropertyLib::Variable::temperature)] = T_i + dT;
+ double const rho_vw1 =
+ property.template value<double>(variable_array, pos, t, dt);
+
+ double const approximated_drho_wv_dT =
+ 0.5 * (rho_vw1 - rho_vw0) / dT;
+
+ double const analytic_drho_wv_dT = property.template dValue<double>(
+ variable_array, MaterialPropertyLib::Variable::temperature, pos,
+ t, dt);
+
+ ASSERT_LE(std::fabs(approximated_drho_wv_dT - analytic_drho_wv_dT),
+ 1e-6)
+ << "for expected derivative of water vapour density with "
+ "respect to temperature "
+ << approximated_drho_wv_dT
+ << " and for computed derivative of water vapour density "
+ "with respect to temperature."
+ << analytic_drho_wv_dT;
+ }
+ }
+
+ // The derivative of the water vapour density with respect of pressure
+ {
+ std::array const pressures = {-1.e+7, -1.e+6, 1.e+5, 1.e+6,
+ 2.e+6, 6.e+6, 1.e+7};
+ std::array const rho_vw_expected = {1.101824, 1.158320, 1.165421,
+ 1.171263, 1.177789, 1.204257,
+ 1.23132};
+
+ for (std::size_t i = 0; i < pressures.size(); ++i)
+ {
+ double const p_i = pressures[i];
+
+ variable_array[static_cast<int>(
+ MaterialPropertyLib::Variable::temperature)] = T;
+ variable_array[static_cast<int>(
+ MaterialPropertyLib::Variable::phase_pressure)] = p_i;
+
+ double const rho_vw =
+ property.template value<double>(variable_array, pos, t, dt);
+
+ ASSERT_LE(std::fabs(rho_vw_expected[i] - rho_vw), 5e-6)
+ << "for expected water vapour density " << rho_vw_expected[i]
+ << " and for computed water vapour density " << rho_vw;
+
+ double const dp = 1.0e-3;
+ variable_array[static_cast<int>(
+ MaterialPropertyLib::Variable::phase_pressure)] = p_i - dp;
+
+ double const rho_vw0 =
+ property.template value<double>(variable_array, pos, t, dt);
+
+ variable_array[static_cast<int>(
+ MaterialPropertyLib::Variable::phase_pressure)] = p_i + dp;
+ double const rho_vw1 =
+ property.template value<double>(variable_array, pos, t, dt);
+
+ double const approximated_drho_wv_dp =
+ 0.5 * (rho_vw1 - rho_vw0) / dp;
+
+ double const analytic_drho_wv_dp = property.template dValue<double>(
+ variable_array, MaterialPropertyLib::Variable::phase_pressure,
+ pos, t, dt);
+
+ ASSERT_LE(std::fabs(approximated_drho_wv_dp - analytic_drho_wv_dp),
+ 1e-7)
+ << "for expected derivative of water vapour density with "
+ "respect to pressure "
+ << approximated_drho_wv_dp
+ << " and for computed derivative of water vapour density "
+ "with respect to pressure."
+ << analytic_drho_wv_dp;
+ }
+ }
+}