diff --git a/bench/fgcs/dolo_fgcs.pqi b/bench/fgcs/dolo_fgcs.pqi
new file mode 100644
index 000000000..39e8553a4
--- /dev/null
+++ b/bench/fgcs/dolo_fgcs.pqi
@@ -0,0 +1,48 @@
+SOLUTION 1
+    units	mol/kgw
+    water	1
+    temperature	25
+    pH	7
+    pe	4
+PURE 1
+    Calcite 0.0 1
+END
+
+RUN_CELLS
+	-cells 1
+END
+
+COPY solution 1 2
+
+#PURE 2
+#    O2g -0.1675 10
+KINETICS 2
+    Calcite
+    -m      0.00207
+    -parms  0.05
+    -tol    1e-10
+    Dolomite
+    -m      0.0
+    -parms  0.01
+    -tol    1e-10
+END
+
+SOLUTION 3
+    units	mol/kgw
+    water 1
+    temp 25 
+    Mg 0.001
+    Cl 0.002
+END
+
+SOLUTION 4
+    units mol/kgw 
+    water 1
+    temp 25
+    Mg 0.002
+    Cl 0.004
+END
+
+RUN_CELLS 
+    -cells 2-4
+END
diff --git a/bench/fgcs/dolo_fgcs_3.R b/bench/fgcs/dolo_fgcs_3.R
new file mode 100644
index 000000000..e8c7a3a4e
--- /dev/null
+++ b/bench/fgcs/dolo_fgcs_3.R
@@ -0,0 +1,116 @@
+rows <- 400
+cols <- 400
+
+grid_def <- matrix(2, nrow = rows, ncol = cols)
+
+# Define grid configuration for POET model
+grid_setup <- list(
+    pqc_in_file = "./dolo_fgcs.pqi",
+    pqc_db_file = "../dolo/phreeqc_kin.dat", # Path to the database file for Phreeqc
+    grid_def = grid_def, # Definition of the grid, containing IDs according to the Phreeqc input script
+    grid_size = c(5, 5), # Size of the grid in meters
+    constant_cells = c() # IDs of cells with constant concentration
+)
+
+bound_def_we <- list(
+    "type" = rep("constant", rows),
+    "sol_id" = rep(1, rows),
+    "cell" = seq(1, rows)
+)
+
+bound_def_ns <- list(
+    "type" = rep("constant", cols),
+    "sol_id" = rep(1, cols),
+    "cell" = seq(1, cols)
+)
+
+diffusion_setup <- list(
+    boundaries = list(
+        "W" = bound_def_we,
+        "E" = bound_def_we,
+        "N" = bound_def_ns,
+        "S" = bound_def_ns
+    ),
+    inner_boundaries = list(
+        "row" = floor(rows / 2),
+        "col" = floor(cols / 2),
+        "sol_id" = c(3)
+    ),
+    alpha_x = 1e-6,
+    alpha_y = 1e-6
+)
+
+check_sign_cal_dol_dht <- function(old, new) {
+    if ((old["Calcite"] == 0) != (new["Calcite"] == 0)) {
+        return(TRUE)
+    }
+    if ((old["Dolomite"] == 0) != (new["Dolomite"] == 0)) {
+        return(TRUE)
+    }
+    return(FALSE)
+}
+
+check_sign_cal_dol_interp <- function(to_interp, data_set) {
+    dht_species <- c(
+        "H" = 3,
+        "O" = 3,
+        "C" = 6,
+        "Ca" = 6,
+        "Cl" = 3,
+        "Mg" = 5,
+        "Calcite" = 4,
+        "Dolomite" = 4
+    )
+    data_set <- as.data.frame(do.call(rbind, data_set), check.names = FALSE, optional = TRUE)
+    names(data_set) <- names(dht_species)
+    cal <- (data_set$Calcite == 0) == (to_interp["Calcite"] == 0)
+    dol <- (data_set$Dolomite == 0) == (to_interp["Dolomite"] == 0)
+
+    cal_dol_same_sig <- cal == dol
+    return(rev(which(!cal_dol_same_sig)))
+}
+
+check_neg_cal_dol <- function(result) {
+    neg_sign <- (result["Calcite"] < 0) || (result["Dolomite"] < 0)
+    return(neg_sign)
+}
+
+# Optional when using Interpolation (example with less key species and custom
+# significant digits)
+
+pht_species <- c(
+    "C"  = 3,
+    "Ca" = 3,
+    "Mg" = 3,
+    "Cl" = 3,
+    "Calcite" = 3,
+    "Dolomite" = 3
+)
+
+
+dht_species <- c(
+        "H" = 3,
+        "O" = 3,
+        "C" = 6,
+        "Ca" = 6,
+        "Cl" = 3,
+        "Mg" = 5,
+        "Calcite" = 4,
+        "Dolomite" = 4)
+
+chemistry_setup <- list(
+    dht_species = dht_species,
+    pht_species = pht_species,
+    hooks = list(
+        dht_fill = check_sign_cal_dol_dht,
+        interp_pre = check_sign_cal_dol_interp,
+        interp_post = check_neg_cal_dol
+    )
+)
+
+# Define a setup list for simulation configuration
+setup <- list(
+    Grid = grid_setup, # Parameters related to the grid structure
+    Diffusion = diffusion_setup, # Parameters related to the diffusion process
+    Chemistry = chemistry_setup # Parameters related to the chemistry process
+)
diff --git a/bench/fgcs/dolo_fgcs_3_rt.R b/bench/fgcs/dolo_fgcs_3_rt.R
new file mode 100644
index 000000000..70731e024
--- /dev/null
+++ b/bench/fgcs/dolo_fgcs_3_rt.R
@@ -0,0 +1,8 @@
+iterations <- 10000
+dt <- 200
+
+list(
+    timesteps = rep(dt, iterations),
+    store_result = TRUE,
+    out_save = out_save = c(1,20,50, seq(100,1000, by=100), seq(2000, iterations, by=1000))
+)