diff --git a/bench/fgcs/barite_fgcs_3.pqi b/bench/fgcs/barite_fgcs_3.pqi
new file mode 100644
index 000000000..a55fdacab
--- /dev/null
+++ b/bench/fgcs/barite_fgcs_3.pqi
@@ -0,0 +1,49 @@
+SOLUTION 1   ## non reactive, at equilibrium with Celestite
+  units	mol/kgw
+  water	1
+  temperature	25
+  pH     7.008 
+  pe    10.798
+  S      6.205e-04
+  Sr     6.205e-04
+END
+
+SOLUTION 2  ## reactive, kinetic Celestite & Barite
+  units	mol/kgw
+  water	1
+  temperature	25
+  pH     7.008 
+  pe    10.798
+  S      6.205e-04
+  Sr     6.205e-04
+KINETICS 2
+  Barite
+    -m 0.00
+    -parms 50.0  # reactive surface area
+    -tol 1e-9
+  Celestite
+    -m 1
+    -parms 20.0  # reactive surface area
+    -tol 1e-9
+END
+
+SOLUTION 3 ## "injection" NW
+  units mol/kgw
+  water 1
+  temperature 25
+  Ba 0.1
+  Cl 0.2
+END
+
+SOLUTION 4 ## "injection" SE
+  units mol/kgw
+  water 1
+  temperature 25
+  Ba 0.2
+  Cl 0.4
+END
+
+
+RUN_CELLS
+  -cells 1 2 3 4
+END
diff --git a/bench/fgcs/barite_fgcs_4.R b/bench/fgcs/barite_fgcs_4.R
new file mode 100644
index 000000000..43902779e
--- /dev/null
+++ b/bench/fgcs/barite_fgcs_4.R
@@ -0,0 +1,109 @@
+## Time-stamp: "Last modified 2024-12-20 15:30:05 delucia"
+
+cols <- 400
+rows <- 400
+
+dim_cols <- 80
+dim_rows <- 80
+
+rad   <- cols/40  ## circle radius (in nodes)
+boxs  <- cols/16  ## box side (in nodes)
+
+cc <- seq(boxs/2, cols - boxs/2, by = boxs)
+length(cc)
+
+colcen <- rep(cc, length(cc))
+rowcen <- rep(cc, each = length(cc))
+centers <- cbind(rowcen, colcen)
+
+mr <- matrix(rep(seq_len(cols), each = rows), byrow = TRUE, nrow = rows)
+mc <- matrix(rep(seq_len(cols), rows), byrow = TRUE, nrow = rows)
+
+tmpl <- lapply(seq_len(nrow(centers)), function(x) which((mr-centers[x, 1])^2 + (mc-centers[x, 2])^2 < rad^2, arr.ind = TRUE))
+
+inds <-  do.call(rbind, tmpl)
+
+grid <- matrix(1, nrow = rows, ncol = cols)
+grid[inds] <- 2
+
+alpha <- matrix( 1e-5, ncol = cols, nrow = rows)
+alpha[inds] <- 1e-8
+
+## PoetUtils::PlotField(grid, cols = cols, rows = rows, contour = FALSE, scale = FALSE, las=1)
+
+## Define grid configuration for POET model
+grid_setup <- list(
+  pqc_in_file    = "./barite_fgcs_3.pqi",
+  pqc_db_file    = "../barite/db_barite.dat",     ## database file
+  grid_def       = grid,                  ## grid definition, IDs according to the Phreeqc input 
+  grid_size      = c(dim_cols, dim_rows), ## grid size in meters
+  constant_cells = c()                    ## IDs of cells with constant concentration
+)
+
+bound_length <- cols/10
+
+bound_N <- list(
+  "type"   = rep("constant", bound_length),
+  "sol_id" = rep(3, bound_length),
+  "cell"   = seq(1, bound_length)
+)
+
+bound_W <- list(
+  "type"   = rep("constant", bound_length),
+  "sol_id" = rep(3, bound_length),
+  "cell"   = seq(1, bound_length)
+)
+bound_E <- list(
+  "type"   = rep("constant", bound_length),
+  "sol_id" = rep(4, bound_length),
+  "cell"   = seq(rows-bound_length+1, rows)
+)
+
+bound_S <- list(
+  "type"   = rep("constant", bound_length),
+  "sol_id" = rep(4, bound_length),
+  "cell"   = seq(cols-bound_length+1, cols)
+)
+
+diffusion_setup <- list(
+  boundaries = list(
+    "W" = bound_W,
+    "N" = bound_N,
+    "E" = bound_E,
+    "S" = bound_S
+  ),
+  alpha_x = alpha,
+  alpha_y = alpha
+)
+
+dht_species <- c(
+  "H"         = 7,
+  "O"         = 7,
+  "Ba"        = 7,
+  "Cl"        = 7,
+  "S"         = 7,
+  "Sr"        = 7,
+  "Barite"    = 4,
+  "Celestite" = 4
+)
+
+pht_species <- c(
+  "Ba"        = 4,
+  "Cl"        = 3,
+  "S"         = 3,
+  "Sr"        = 3,
+  "Barite"    = 2,
+  "Celestite" = 2
+)
+
+chemistry_setup <- list(
+  dht_species = dht_species,
+  pht_species = pht_species
+)
+
+## 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
+)
diff --git a/bench/fgcs/barite_fgcs_4_rt.R b/bench/fgcs/barite_fgcs_4_rt.R
new file mode 100644
index 000000000..601556c56
--- /dev/null
+++ b/bench/fgcs/barite_fgcs_4_rt.R
@@ -0,0 +1,10 @@
+iterations <- 5000
+
+dt <- 2000
+
+list(
+    timesteps = rep(dt, iterations),
+    store_result = TRUE,
+    out_save = c(seq(1, 10), seq(10, 100, by= 10), seq(200, iterations, by=100))
+    ## out_save = seq_len(iterations)
+)