diff --git a/R_lib/init_r_lib.R b/R_lib/init_r_lib.R
index 672e372cc..c9b096057 100644
--- a/R_lib/init_r_lib.R
+++ b/R_lib/init_r_lib.R
@@ -28,3 +28,14 @@ pqc_to_grid <- function(pqc_in, grid) {
 
     return(res_df)
 }
+
+resolvePqcBound <- function(pqc_mat, transport_spec, id) {
+    df <- as.data.frame(pqc_mat, check.names = FALSE)
+    value <- df[df$ID == id, transport_spec]
+
+    if (is.nan(value)) {
+        value <- 0
+    }
+
+    return(value)
+}
diff --git a/src/Init/DiffusionInit.cpp b/src/Init/DiffusionInit.cpp
index a30b66ec6..93046575d 100644
--- a/src/Init/DiffusionInit.cpp
+++ b/src/Init/DiffusionInit.cpp
@@ -26,49 +26,6 @@ const std::map<std::uint8_t, std::string> tug_side_mapping = {
     {tug::BC_SIDE_TOP, "N"},
     {tug::BC_SIDE_BOTTOM, "S"}};
 
-static Rcpp::List parseBoundaries2D(const Rcpp::List &boundaries_list,
-                                    std::uint32_t n_cols,
-                                    std::uint32_t n_rows) {
-  Rcpp::List out_list;
-
-  for (const auto &side : tug_side_mapping) {
-    Rcpp::NumericVector type =
-        (side.first == tug::BC_SIDE_RIGHT || side.first == tug::BC_SIDE_LEFT)
-            ? Rcpp::NumericVector(n_rows, tug::BC_TYPE_CLOSED)
-            : Rcpp::NumericVector(n_cols, tug::BC_TYPE_CLOSED);
-
-    Rcpp::NumericVector value =
-        (side.first == tug::BC_SIDE_RIGHT || side.first == tug::BC_SIDE_LEFT)
-            ? Rcpp::NumericVector(n_rows, 0)
-            : Rcpp::NumericVector(n_cols, 0);
-
-    if (boundaries_list.containsElementNamed(side.second.c_str())) {
-      const Rcpp::List mapping = boundaries_list[side.second];
-
-      const Rcpp::NumericVector cells = mapping["cells"];
-      const Rcpp::NumericVector values = mapping["sol_id"];
-
-      if (cells.size() != values.size()) {
-        throw std::runtime_error("Boundary [" + side.second +
-                                 "] cells and values are not the same length");
-      }
-    }
-
-    out_list[side.second] = Rcpp::List::create(Rcpp::Named("type") = type,
-                                               Rcpp::Named("value") = value);
-  }
-}
-
-static inline SEXP_TYPE get_datatype(const SEXP &input) {
-  Rcpp::Function check_list("is.list");
-
-  if (Rcpp::as<bool>(check_list(input))) {
-    return SEXP_IS_LIST;
-  } else {
-    return SEXP_IS_VEC;
-  }
-}
-
 static std::vector<TugType> colMajToRowMaj(const Rcpp::NumericVector &vec,
                                            std::uint32_t n_cols,
                                            std::uint32_t n_rows) {
@@ -90,6 +47,70 @@ static std::vector<TugType> colMajToRowMaj(const Rcpp::NumericVector &vec,
   }
 }
 
+Rcpp::List InitialList::resolveBoundaries(const Rcpp::List &boundaries_list) {
+  Rcpp::List bound_list;
+  Rcpp::Function resolve_R("resolvePqcBound");
+
+  for (const auto &species : this->transport_names) {
+    Rcpp::List spec_list;
+
+    for (const auto &side : tug_side_mapping) {
+      Rcpp::NumericVector c_type =
+          (side.first == tug::BC_SIDE_RIGHT || side.first == tug::BC_SIDE_LEFT)
+              ? Rcpp::NumericVector(n_rows, tug::BC_TYPE_CLOSED)
+              : Rcpp::NumericVector(n_cols, tug::BC_TYPE_CLOSED);
+
+      Rcpp::NumericVector c_value =
+          (side.first == tug::BC_SIDE_RIGHT || side.first == tug::BC_SIDE_LEFT)
+              ? Rcpp::NumericVector(n_rows, 0)
+              : Rcpp::NumericVector(n_cols, 0);
+
+      if (boundaries_list.containsElementNamed(side.second.c_str())) {
+        const Rcpp::List mapping = boundaries_list[side.second];
+
+        const Rcpp::NumericVector cells = mapping["cell"];
+        const Rcpp::NumericVector values = mapping["sol_id"];
+        const Rcpp::CharacterVector type_str = mapping["type"];
+
+        if (cells.size() != values.size()) {
+          throw std::runtime_error(
+              "Boundary [" + side.second +
+              "] cells and values are not the same length");
+        }
+
+        for (std::size_t i = 0; i < cells.size(); i++) {
+          if (type_str[i] == "closed") {
+            c_type[cells[i] - 1] = tug::BC_TYPE_CLOSED;
+          } else if (type_str[i] == "constant") {
+            c_type[cells[i] - 1] = tug::BC_TYPE_CONSTANT;
+            c_value[cells[i] - 1] = Rcpp::as<TugType>(
+                resolve_R(this->phreeqc_mat, Rcpp::wrap(species), values[i]));
+          } else {
+            throw std::runtime_error("Unknown boundary type");
+          }
+        }
+      }
+
+      spec_list[side.second] = Rcpp::List::create(
+          Rcpp::Named("type") = c_type, Rcpp::Named("value") = c_value);
+    }
+
+    bound_list[species] = spec_list;
+  }
+
+  return bound_list;
+}
+
+static inline SEXP_TYPE get_datatype(const SEXP &input) {
+  Rcpp::Function check_list("is.list");
+
+  if (Rcpp::as<bool>(check_list(input))) {
+    return SEXP_IS_LIST;
+  } else {
+    return SEXP_IS_VEC;
+  }
+}
+
 static Rcpp::List parseAlphas(const SEXP &input,
                               const std::vector<std::string> &transport_names,
                               std::uint32_t n_cols, std::uint32_t n_rows) {
@@ -132,23 +153,25 @@ static Rcpp::List parseAlphas(const SEXP &input,
   return out_list;
 }
 void InitialList::initDiffusion(const Rcpp::List &diffusion_input) {
-  // const Rcpp::List &boundaries =
-  //     diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::BOUNDARIES)];
+  const Rcpp::List &boundaries =
+      diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::BOUNDARIES)];
   const SEXP &alpha_x =
       diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::ALPHA_X)];
   const SEXP &alpha_y =
       diffusion_input[DIFFU_MEMBER_STR(DiffusionMembers::ALPHA_Y)];
 
+  this->boundaries = resolveBoundaries(boundaries);
+
   this->alpha_x =
       parseAlphas(alpha_x, this->transport_names, this->n_cols, this->n_rows);
 
   this->alpha_y =
       parseAlphas(alpha_y, this->transport_names, this->n_cols, this->n_rows);
 
-  R["alpha_x"] = this->alpha_x;
-  R["alpha_y"] = this->alpha_y;
+  // R["alpha_x"] = this->alpha_x;
+  // R["alpha_y"] = this->alpha_y;
 
-  R.parseEval("print(alpha_x)");
-  R.parseEval("print(alpha_y)");
+  // R.parseEval("print(alpha_x)");
+  // R.parseEval("print(alpha_y)");
 }
 } // namespace poet
\ No newline at end of file
diff --git a/src/Init/InitialList.hpp b/src/Init/InitialList.hpp
index 3c37f84ee..fe1afcc59 100644
--- a/src/Init/InitialList.hpp
+++ b/src/Init/InitialList.hpp
@@ -93,6 +93,8 @@ private:
   }
 
   void initDiffusion(const Rcpp::List &diffusion_input);
+  Rcpp::List resolveBoundaries(const Rcpp::List &boundaries_list);
+
   Rcpp::List boundaries;
 
   Rcpp::List alpha_x;