diff --git a/ChemistryLib/PhreeqcIOData/AqueousSolution.cpp b/ChemistryLib/PhreeqcIOData/AqueousSolution.cpp
index fe0a5d0f18c6efa705f14093c12e7effba782c08..8fdd0635e3dcb696e5576ff1668f07ac073a82a5 100644
--- a/ChemistryLib/PhreeqcIOData/AqueousSolution.cpp
+++ b/ChemistryLib/PhreeqcIOData/AqueousSolution.cpp
@@ -12,6 +12,7 @@
 
 #include "AqueousSolution.h"
 #include "ChemistryLib/Common/ChargeBalance.h"
+#include "MeshLib/PropertyVector.h"
 
 namespace ChemistryLib
 {
diff --git a/ChemistryLib/PhreeqcIOData/AqueousSolution.h b/ChemistryLib/PhreeqcIOData/AqueousSolution.h
index 43532f62d19ce6a948c87effb457f0b72a6c916f..000c1e76667622770f6568c84cf5fe1dd8d414b6 100644
--- a/ChemistryLib/PhreeqcIOData/AqueousSolution.h
+++ b/ChemistryLib/PhreeqcIOData/AqueousSolution.h
@@ -19,6 +19,12 @@
 #include "MathLib/LinAlg/MatrixVectorTraits.h"
 #include "Output.h"
 
+namespace MeshLib
+{
+template <typename PROP_VAL_TYPE>
+class PropertyVector;
+}
+
 namespace ChemistryLib
 {
 enum class ChargeBalance;
@@ -45,6 +51,7 @@ struct AqueousSolution
     AqueousSolution(double temperature_,
                     double pressure_,
                     double pe_,
+                    MeshLib::PropertyVector<double>* pe_,
                     std::vector<Component>&& components_,
                     ChargeBalance charge_balance_,
                     std::size_t const num_chemical_systems_)
@@ -63,8 +70,8 @@ struct AqueousSolution
 
     double temperature;
     double pressure;
-    double pe;
     std::unique_ptr<GlobalVector> pH;
+    MeshLib::PropertyVector<double>* pe;
     std::vector<Component> components;
     ChargeBalance const charge_balance;
 };
diff --git a/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.cpp b/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.cpp
index 70a7e2c6157b7689adb68b372e72d35912731303..91ba7f9015768c1bec1fde5aea1ecebea13fbb72 100644
--- a/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.cpp
+++ b/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.cpp
@@ -13,12 +13,16 @@
 #include "BaseLib/ConfigTree.h"
 #include "ChemistryLib/Common/CreateChargeBalance.h"
 #include "CreateSolutionComponent.h"
+#include "MeshLib/Mesh.h"
 
 namespace ChemistryLib
 {
 namespace PhreeqcIOData
 {
-AqueousSolution createAqueousSolution(BaseLib::ConfigTree const& config)
+std::unique_ptr<AqueousSolution> createAqueousSolution(
+    BaseLib::ConfigTree const& config,
+    MeshLib::Mesh const& mesh,
+    MeshLib::PropertyVector<std::size_t> const& chemical_system_map)
 {
     //! \ogs_file_param{prj__chemical_system__solution__temperature}
     auto const temperature = config.getConfigParameter<double>("temperature");
@@ -27,13 +31,31 @@ AqueousSolution createAqueousSolution(BaseLib::ConfigTree const& config)
     auto const pressure = config.getConfigParameter<double>("pressure");
 
     //! \ogs_file_param{prj__chemical_system__solution__pe}
-    auto const pe = config.getConfigParameter<double>("pe");
+    auto const pe0 = config.getConfigParameter<double>("pe");
 
-    auto components = createSolutionComponents(config);
+    auto pe = MeshLib::getOrCreateMeshProperty<double>(
+        const_cast<MeshLib::Mesh&>(mesh), "pe", MeshLib::MeshItemType::Node, 1);
+
+    std::fill(std::begin(*pe),
+              std::end(*pe),
+              std::numeric_limits<double>::quiet_NaN());
+
+    std::for_each(
+        chemical_system_map.begin(),
+        chemical_system_map.end(),
+        [&pe, pe0](auto const& global_id) { (*pe)[global_id] = pe0; });
+
+    auto components =
+        createSolutionComponents(config, mesh.getNumberOfBaseNodes());
 
     auto charge_balance = createChargeBalance(config);
 
-    return {temperature, pressure, pe, std::move(components), charge_balance};
+    return std::make_unique<AqueousSolution>(temperature,
+                                             pressure,
+                                             pe,
+                                             std::move(components),
+                                             charge_balance,
+                                             mesh.getNumberOfBaseNodes());
 }
 }  // namespace PhreeqcIOData
 }  // namespace ChemistryLib
diff --git a/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.h b/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.h
index 8e0d680f276eca3a44bf227a0a6c7c9576db7a21..d83f595b92433cc704bbb7390dc76250b60de46e 100644
--- a/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.h
+++ b/ChemistryLib/PhreeqcIOData/CreateAqueousSolution.h
@@ -10,17 +10,30 @@
 
 #pragma once
 
+#include <memory>
+
 namespace BaseLib
 {
 class ConfigTree;
 }
 
+namespace MeshLib
+{
+class Mesh;
+
+template <typename PROP_VAL_TYPE>
+class PropertyVector;
+}  // namespace MeshLib
+
 namespace ChemistryLib
 {
 namespace PhreeqcIOData
 {
 struct AqueousSolution;
 
-AqueousSolution createAqueousSolution(BaseLib::ConfigTree const& config);
+std::unique_ptr<AqueousSolution> createAqueousSolution(
+    BaseLib::ConfigTree const& config,
+    MeshLib::Mesh const& mesh,
+    MeshLib::PropertyVector<std::size_t> const& chemical_system_map);
 }  // namespace PhreeqcIOData
 }  // namespace ChemistryLib