diff --git a/bench/surfex/CMakeLists.txt b/bench/surfex/CMakeLists.txt
index 8b3f44d5b..c5a605c6c 100644
--- a/bench/surfex/CMakeLists.txt
+++ b/bench/surfex/CMakeLists.txt
@@ -1,8 +1,8 @@
 install(FILES
     ExBase.pqi
     ex.R
-    #surfex.R
-    #SurfExBase.pqi
+    surfex.R
+    SurfExBase.pqi
     SMILE_2021_11_01_TH.dat
 DESTINATION
     share/poet/bench/surfex
diff --git a/bench/surfex/surfex.R b/bench/surfex/surfex.R
index c2bd714c8..37a5e0429 100644
--- a/bench/surfex/surfex.R
+++ b/bench/surfex/surfex.R
@@ -1,4 +1,4 @@
-## Time-stamp: "Last modified 2023-04-13 17:11:52 mluebke"
+## Time-stamp: "Last modified 2023-04-17 15:48:21 mluebke"
 
 database <- normalizePath("../share/poet/bench/surfex/SMILE_2021_11_01_TH.dat")
 input_script <- normalizePath("../share/poet/bench/surfex/SurfExBase.pqi")
@@ -103,8 +103,11 @@ boundary <- list(
 
 diffu_list <- names(alpha_diffu)
 
+vecinj <- do.call(rbind.data.frame, vecinj_diffu)
+names(vecinj) <- names(init_cell)
+
 diffusion <- list(
-  init = init_cell,
+  init = as.data.frame(init_cell, check.names = FALSE),
   vecinj = vecinj,
 #  vecinj_inner = vecinj_inner,
   vecinj_index = boundary,
diff --git a/ext/phreeqcrm b/ext/phreeqcrm
index ade4f5660..6a4c14c28 160000
--- a/ext/phreeqcrm
+++ b/ext/phreeqcrm
@@ -1 +1 @@
-Subproject commit ade4f56605b8d0e81b7ccad510f2db21f46b5ff9
+Subproject commit 6a4c14c2857f0ad3558af3f1e1aa7555a805e8e2
diff --git a/include/poet/ChemistryModule.hpp b/include/poet/ChemistryModule.hpp
index f37f1f2de..fc096c862 100644
--- a/include/poet/ChemistryModule.hpp
+++ b/include/poet/ChemistryModule.hpp
@@ -385,9 +385,9 @@ protected:
   uint32_t prop_count = 0;
   std::vector<std::string> prop_names;
   std::vector<double> field;
-  static constexpr uint32_t MODULE_COUNT = 4;
+  static constexpr int MODULE_COUNT = 5;
 
-  std::array<std::uint32_t, 4> speciesPerModule{};
+  std::array<std::uint32_t, MODULE_COUNT> speciesPerModule{};
 };
 } // namespace poet
 
diff --git a/src/ChemistryModule/ChemistryModule.cpp b/src/ChemistryModule/ChemistryModule.cpp
index d46ddf60a..3357e0397 100644
--- a/src/ChemistryModule/ChemistryModule.cpp
+++ b/src/ChemistryModule/ChemistryModule.cpp
@@ -68,6 +68,7 @@ void poet::ChemistryModule::RunInitFile(const std::string &input_script_path) {
   char exchange = (speciesPerModule[1] == 0 ? -1 : 1);
   char kinetics = (speciesPerModule[2] == 0 ? -1 : 1);
   char equilibrium = (speciesPerModule[3] == 0 ? -1 : 1);
+  char surface = (speciesPerModule[4] == 0 ? -1 : 1);
 
 #ifdef POET_USE_PRM
   std::vector<int> ic1;
@@ -77,7 +78,7 @@ void poet::ChemistryModule::RunInitFile(const std::string &input_script_path) {
     ic1[i] = 1;                   // Solution 1
     ic1[nxyz + i] = equilibrium;  // Equilibrium 1
     ic1[2 * nxyz + i] = exchange; // Exchange none
-    ic1[3 * nxyz + i] = -1;       // Surface none
+    ic1[3 * nxyz + i] = surface;  // Surface none
     ic1[4 * nxyz + i] = -1;       // Gas phase none
     ic1[5 * nxyz + i] = -1;       // Solid solutions none
     ic1[6 * nxyz + i] = kinetics; // Kinetics 1
@@ -92,7 +93,7 @@ void poet::ChemistryModule::RunInitFile(const std::string &input_script_path) {
     ic1[i] = 1;                   // Solution 1
     ic1[nxyz + i] = equilibrium;  // Equilibrium 1
     ic1[2 * nxyz + i] = exchange; // Exchange none
-    ic1[3 * nxyz + i] = -1;       // Surface none
+    ic1[3 * nxyz + i] = surface;  // Surface none
     ic1[4 * nxyz + i] = -1;       // Gas phase none
     ic1[5 * nxyz + i] = -1;       // Solid solutions none
     ic1[6 * nxyz + i] = kinetics; // Kinetics 1
@@ -104,7 +105,7 @@ void poet::ChemistryModule::RunInitFile(const std::string &input_script_path) {
 
 #ifndef POET_USE_PRM
 void poet::ChemistryModule::mergeFieldWithModule(const SingleCMap &input_map,
-                                               std::uint32_t n_cells) {
+                                                 std::uint32_t n_cells) {
 
   if (is_master) {
     int f_type = CHEM_INIT_SPECIES;
@@ -171,16 +172,8 @@ void poet::ChemistryModule::mergeFieldWithModule(const SingleCMap &input_map,
 
     const std::vector<std::string> ess_names = old_prop_names;
 
-    for (const auto &name : this->prop_names) {
-      auto it_find = std::find(ess_names.begin(), ess_names.end(), name);
-      double value;
-      if (it_find != ess_names.end()) {
-        int index = it_find - ess_names.begin();
-        value = init_values[index];
-      } else {
-        auto map_it = input_map.find(name);
-        value = map_it->second;
-      }
+    for (int i = 0; i < prop_names.size(); i++) {
+      double value = init_values[i];
       std::vector<double> curr_vec(n_cells, value);
       this->field.insert(field.end(), curr_vec.begin(), curr_vec.end());
     }