diff --git a/hydrus/README.md b/hydrus/README.md
new file mode 100644
index 0000000000000000000000000000000000000000..24b7691cb69efdd8c2b8cba036aea78bd9b2796e
--- /dev/null
+++ b/hydrus/README.md
@@ -0,0 +1,2 @@
+# TW modeling toolkit
+
diff --git a/hydrus/example1.zip b/hydrus/example1.zip
new file mode 100644
index 0000000000000000000000000000000000000000..5000cec0df6c56020d513bead5faa1898404e5d7
Binary files /dev/null and b/hydrus/example1.zip differ
diff --git a/hydrus/example2.zip b/hydrus/example2.zip
new file mode 100644
index 0000000000000000000000000000000000000000..36b78da432c2e7dba1d45a37802f2eaac26ac80c
Binary files /dev/null and b/hydrus/example2.zip differ
diff --git a/hydrus/file_manip.R b/hydrus/file_manip.R
new file mode 100644
index 0000000000000000000000000000000000000000..fd195fa07671c78274679606b14136e0a164e518
--- /dev/null
+++ b/hydrus/file_manip.R
@@ -0,0 +1,188 @@
+# Nicolas Forquet, 2015
+
+new.hydrus.subfolder = function(new.name,reference.folder,hydrus.folder){
+  # new.name : name of the new hydrus subfolder
+  # reference.folder : folder containing all new hydrus subfolder
+  # hydrus.folder : folder containing the reference model obtained with the GUI
+  # NOTE: the reference folder must contain a copy of the Hydrus batch executable file : H2D_Calc.exe
+  
+  # get the content list from the reference folder
+  files2copy <- list.files(path = hydrus.folder)
+  # create a new folder into the reference folder
+  setwd(reference.folder)
+  dir.create(new.name)
+  # copy all input files
+  file.copy("H2D_Calc.exe",paste(new.name,"/H2D_Calc.exe",sep=""))
+  for (i in 1:length(files2copy)){
+    file.copy(paste(hydrus.folder,files2copy[i],sep = "/"),paste(new.name,files2copy[i],sep = "/"))
+  }
+}
+
+dummy.white.space = function(char.line){
+  temp = unlist(strsplit(char.line," ",fixed = TRUE))
+  test.temp = (temp =="")
+  result.vector = temp[!test.temp]
+  return(result.vector)
+}
+
+read.meshtria = function(reference.folder){
+  setwd(reference.folder)
+  conn = file("meshtria.txt",open="r")
+  linn = readLines(conn)
+  close(conn)
+  
+  
+  nb.nodes = as.numeric(dummy.white.space(linn[1])[2])
+  nb.mesh = as.numeric(dummy.white.space(linn[1])[4])
+  coord.nodes = data.frame(x = rep(0,nb.nodes), z = rep(0,nb.nodes))
+  edges = data.frame(i = rep(0,nb.mesh), j = rep(0,nb.mesh), k = rep(0, nb.mesh))
+  
+  for (i in 1:nb.nodes){
+    coord.nodes$x[i] = as.numeric(dummy.white.space(linn[i+1])[2])
+    coord.nodes$y[i] = as.numeric(dummy.white.space(linn[i+1])[3])
+  }
+  end.of.table = i+1
+  for (m in 1:nb.mesh){
+    edges$i[m] = as.numeric(dummy.white.space(linn[m+end.of.table+3])[2])
+    edges$j[m] = as.numeric(dummy.white.space(linn[m+end.of.table+3])[3])
+    edges$k[m] = as.numeric(dummy.white.space(linn[m+end.of.table+3])[4])
+  }
+  
+  output = list(nb.nodes = nb.nodes, nb.mesh = nb.mesh, coord.nodes = coord.nodes, edges = edges, linn = linn)
+  return(output)
+}
+
+write_domain = function(reference.folder, file.name.out = "empty", Axz = "empty", Bxz = "empty", nb.col = 11, header.size = 6){
+  # set working directory to reference folder
+  setwd(reference.folder)
+  # Hydrus has been developed with Windows. This is bad. What is even worse is that Windows is not case sensitive and therefore Simunek sometimes names its file in upper case, sometimes in lower case. Therefore I introduced a step that list the files having the name we are looking for independly of the case
+  file.name.in <- list.files(pattern = "domain.dat", ignore.case = TRUE)
+  
+  # Then we open the file and extract the lines content and put in linn.
+  conn <- file(file.name.in,open="r")
+  linn <- readLines(conn)
+  close(conn) 
+  
+  # test on model arguments
+  if ((is.character(Axz) == TRUE) & (is.character(Bxz) == TRUE)){
+    stop("At least Axz or Bxz must be given")
+  }
+  
+  Axz.log <- is.character(Axz)
+  Bxz.log <- is.character(Bxz)
+  if (Axz.log == FALSE){
+    nb.nodes <- length(Axz)
+  } else {
+    nb.nodes <- length(Bxz)  
+  }
+   
+  # prealocation of matrix that contains data from the domain.dat file
+  node.mat = matrix(rep(0,nb.col*nb.nodes),ncol=nb.col)
+  
+  # filling up node.map line by line starting after the header.size
+  for (i in 1:nb.nodes){
+    node.mat[i,] = as.numeric(dummy.white.space(linn[header.size+i]))
+  }
+  
+  if (Axz.log == FALSE){
+    node.mat[,7] = Axz
+  }
+  if (Bxz.log == FALSE){
+    node.mat[,8] = Bxz  
+  }
+  
+  new_linn = linn
+  for (i in 1:nb.nodes){
+    new_linn[i+header.size] = paste(as.character(node.mat[i,]),collapse=" ")
+  }
+  
+  new_linn[(header.size+nb.nodes+1):length(linn)] = linn[(header.size+nb.nodes+1):length(linn)]
+  
+  if (file.name.out == "empty"){
+    file.name.out <- file.name.in
+  }
+  
+  conn = file(file.name.out,open="w")
+  writeLines(new_linn,conn)
+  close(conn) 
+  
+  return(new_linn)
+}
+
+level_01.creator = function(new.name,reference.folder){
+  # new.name : name of the new hydrus subfolder
+  # reference.folder : folder containing all new hydrus subfolder
+  level01.file = file(paste(new.name,"/Level_01.dir",sep=""),"w")
+  working.directory = paste(reference.folder,"/",new.name,sep="")
+  writeLines(working.directory,level01.file)
+  close(level01.file)
+}
+
+modify.selector = function(new.name,hydrus.folder,data,header.size = 24, end.size = 65){
+  # nombre de materiaux dans le modele
+  nb.mat <-dim(data)[1]
+  #read old file to get data written before and after the material matrix
+  selector.file = file(paste(hydrus.folder,"/SELECTOR.IN",sep=""),"r")
+  #store the header.size first lines (no transformation will be done on those lines)
+  before.text = readLines(selector.file,n=header.size)
+  # read data with material properties (data not save)
+  readLines(selector.file,n=dim(data)[1])
+  # store file content after material properties
+  after.text = readLines(selector.file,n=end.size)
+  #test.data = list(before = before.text,mat1 = mat1.text,mid = mid.text, mat2 = mat2.text, after = after.text)
+  close(selector.file)
+  
+  #create the new selector.in in the appropriate folder
+  selector.new.file = file(paste(new.name,"/SELECTOR.IN",sep=""),"w")
+  # write text block before material properties
+  writeLines(before.text,selector.new.file)
+  # enter the new material properties
+  write.table(data,selector.new.file,row.names=FALSE,col.names=FALSE,sep="\t")
+  # write text block after material properties
+  writeLines(after.text,selector.new.file)
+  # close file connection
+  close(selector.new.file)
+}
+
+modify.hydraulic.parameters = function(new.name,hydrus.folder,data,hydrus.version="1D"){
+  # so far this function as only be tested for hydrus 1D
+  #find the header line number corresponding to hydraulic parameters
+  selector.file <- file(paste(hydrus.folder,"/SELECTOR.IN",sep=""),"r")
+  file.content <- readLines(selector.file,n=-1)
+  nb.line.start <- which(grepl("thr",file.content))
+  close(selector.file)
+  # number of model materials
+  nb.mat <-dim(data)[1]
+  #store the lines before the hydraulic parameters block (no transformation will be done on those lines)
+  before.text = file.content[1:nb.line.start]
+  #store the lines after the hydraulic parameters block (no transformation will be done on those lines)
+  after.text = file.content[nb.line.start+nb.mat+1:length(file.content)]
+  #create the new selector.in in the appropriate folder
+  selector.new.file = file(paste(new.name,"/SELECTOR.IN",sep=""),"w")
+  # write text block before material properties
+  writeLines(before.text,selector.new.file)
+  # enter the new material properties
+  write.table(data,selector.new.file,row.names=FALSE,col.names=FALSE,sep="\t")
+  # write text block after material properties
+  writeLines(after.text,selector.new.file)
+  # close file connection
+  close(selector.new.file)
+}
+
+read.binary = function(filename,nb.elts,nb.ts){
+  fileSize <- file.info(filename)$size
+  raw <- readBin(filename,what="raw",n=fileSize)
+  if (length(raw) %% 4 !=0){
+    stop('File format assumption error')
+  }
+  nbrOfRecords <- length(raw)/4
+  unsort.data = readBin(con=raw,what="double",size = 4,n = nbrOfRecords,endian="little")
+  time.ts = rep(0,nb.ts)
+  data.ts = matrix(rep(0,(nb.ts*nb.elts)),ncol = nb.ts)
+  for(i in 1:nb.ts){
+    time.ts[i] = unsort.data[i+nb.elts*(i-1)]
+    data.ts[1:nb.elts,i] = unsort.data[(nb.elts*(i-1)+i+1):(i*(nb.elts+1))]
+  }
+  output = list(time.ts = time.ts, data.ts = data.ts)
+  return(output)
+}
\ No newline at end of file
diff --git a/hydrus/loop_example.Rmd b/hydrus/loop_example.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..05995fac4583bbd22504af96142cc12613412050
--- /dev/null
+++ b/hydrus/loop_example.Rmd
@@ -0,0 +1,76 @@
+---
+title: "How to run several hydrus models in batch mode"
+output: html_notebook
+---
+
+# Prerequires
+
+1. You must have created a reference model using the GUI. Create a **permanent model** within a folder located in your User directory.
+2. Make sure to uncheck the box *Hit Enter at the End?* located in the *Print Information* window.
+3. Copy the appropriate exe file from *Program Files* directory (Hydrus is usually installed within a folder named *PC-progress*) to your model folder.
+
+For example, if you want to run a dual porosity 1D model, you must copy the file named *H1D_Dual.exe*.
+
+# load functions and define paths
+
+You need to specify the location of the rhydrus folder:
+```{r}
+rhydrus.folder <- "C:/Users/forquet/Documents/script_R/rhydrus"
+```
+
+To load all necessary functions, source these two files:
+```{r}
+filename <- paste0(rhydrus.folder,"/file_manip.R")
+source(filename)
+filename <- paste0(rhydrus.folder,"/simulation_fcn.R")
+source(filename)
+```
+
+The main directory, i.e. where you saved your model and where a new folder will be created for each run you will set up.
+```{r}
+main.directory <- "C:/Users/forquet/Documents/script_R/rhydrus_demo/PF_1pulse"
+setwd(main.directory)
+```
+
+Indicate the number of run to be carried out:
+```{r}
+nb.runs <- 2
+```
+
+As explained earlier, each simulation results will be placed in a new folder. It's name will be composed of a common root and the number of the run. The common root is named *generic.name*.
+```{r}
+generic.name <- "PF_Model_"
+```
+
+With the main directory, you must have placed a gui file ending with an extension *.h1d* and a folder containing all input and output files. The name of this folder must be saved within the variable *hydrus.folder*.
+```{r}
+hydrus.folder <- "DPo_1pulse_modifpdt_180108"
+```
+
+Now you must specify the new hydraulic parameters you want to use for each run. in this example we run a dual porosity model with 2 materials. Therefore you must create a matrix for each material with the number of line equal to the number of run and the number of columns equal to the number of parameters (9 in dual porosity model).
+```{r}
+# material 1
+deposit <- matrix(c(0, 0.05, 0.12, 1.8, 5, 0.5, 0.64, 0.79, 0.05, 0, 0.05, 0.12, 1.8, 5, 0.5, 0.64, 0.79, 0.05),ncol=9,byrow=TRUE)
+
+# material 2
+gravel <- matrix(c(0, 0.1, 0.572, 2.917, 300, 0.5, 0.052, 0.4, 0, 0, 0.1, 0.572, 2.917, 300, 0.5, 0.052, 0.4, 0),ncol=9,byrow=TRUE)
+```
+
+A list is created made of a data.frame for each run
+```{r}
+new.data <- list()
+for (i in 1:nb.runs){
+  new.data[[i]] <- rbind(deposit[i,],gravel[i,])
+}
+```
+
+Next we use the *simulation.setup* function to create all the folders and to copy input files.
+```{r}
+simulation.setup(nb.runs,generic.name,main.directory,hydrus.folder,new.data,hydrus.version="1D")
+```
+
+Finally, we run the simulations using the function *simulation.run*. This function output the computing time:
+```{r eval=FALSE, include=FALSE}
+output <- simulation.run(2,paste0(rhydrus.folder,"/processor_load.bat"),main_directory,generic.name,hydrus.calc="H1D_Dual.exe")
+```
+
diff --git a/hydrus/main_script.R b/hydrus/main_script.R
new file mode 100644
index 0000000000000000000000000000000000000000..d37eea0cfa4b83b44251661a3091a9ae55a74a08
--- /dev/null
+++ b/hydrus/main_script.R
@@ -0,0 +1,76 @@
+# specify the location of the rhydrus folder:
+rhydrus.folder <- "C:/Users/forquet/Documents/script_R/rhydrus"
+# source the file with functions dedicated to file manipulation
+filename <- paste0(rhydrus.folder,"/file_manip.R")
+source(filename)
+
+# This function process the output of the script .bat that estimates processor usage.
+output.bat <- function(output){
+  usage <- rep(0,5)
+  for (i in 1:5){
+    usage[i] <- as.numeric(strsplit(output[4+i],split='"',fixed=TRUE)[[1]][4])    
+  }
+  mean_usage <- mean(usage)
+  return(mean_usage)
+}
+
+# repository that contains all Hydrus folders
+main.directory <- "C:/Users/forquet/Documents/script_R/rhydrus_demo/PF_1pulse"
+
+# set the working directory to main.directory 
+setwd(main.directory)
+
+# nb of simulation to run
+nb.runs <- 1
+
+# generic name = the name that will be given to the model that will be run iteratively
+# hydrus.folder = the original hydrus folder generated using the GUI
+# VERY IMPORTANT = the reference folder MUST contain the HYDRUS exe file !!!
+# The option press enter at the end must have been unchecked in the GUI
+generic.name <- "PF_Model_"
+hydrus.folder <- "DPo_1pulse_modifpdt_180108"
+
+## van Genuchten parameters (here you can enter matrices for the different cases you want to test)
+# gravel parameters
+deposit = matrix(c(0, 0.05, 0.12, 1.8, 5, 0.5, 0.64, 0.79, 0.05),ncol=9)
+
+# sand parameters
+gravel = matrix(c(0, 0.1, 0.572, 2.917, 300, 0.5, 0.052, 0.4, 0),ncol=9)
+
+# create all simulation folders, copy files and modify the one containing the van Genuchten parameters
+for (i in 1:nb.runs){
+  # current simulation name
+  new.name <- paste(generic.name,as.character(i),sep = "")
+  # create a copy of the reference model in a folder named as the current simulation
+  new.hydrus.subfolder(new.name,main.directory,hydrus.folder)
+  # create the level_01 file
+  level_01.creator(new.name,main.directory)
+  # modification of the selector.in file
+  new.data <- rbind(deposit,gravel)
+  modify.hydraulic.parameters(new.name,hydrus.folder,new.data,hydrus.version="1D")
+}
+
+# computation loop
+start.time <- Sys.time()
+sim.counter <- 0
+
+while (sim.counter < nb.runs){
+  # how busy the processor is?
+  temp <- system('C:/Users/forquet/Documents/script_R/rhydrus/processor_load.bat',intern = TRUE)
+  processor.usage <- output.bat(temp)
+  if (processor.usage < 75){
+    # if processor usage is less than 75%, start a new run
+    sim.counter <- sim.counter + 1
+    # go to the simulation working directory
+    setwd(paste(main.directory,"/",generic.name,as.character(sim.counter),sep=""))
+    # run Hydrus
+    system("H1D_Dual.exe", wait = FALSE)
+    # print the current simulation number
+    print(sim.counter)
+    # go back to main directory
+    setwd(main.directory)
+  }
+  # wait a minute before restarting the loop
+  Sys.sleep(60)
+}
+end.time <- Sys.time()
\ No newline at end of file
diff --git a/hydrus/postprocessing_fcn.R b/hydrus/postprocessing_fcn.R
new file mode 100644
index 0000000000000000000000000000000000000000..86a4f518f66ea8c037d616c26c259498173efa22
--- /dev/null
+++ b/hydrus/postprocessing_fcn.R
@@ -0,0 +1,78 @@
+library(caTools)
+
+read.solute.out <- function(filename){
+  # lecture du fichier ligne par ligne
+  conn <- file(filename,open="r")
+  linn = readLines(conn)
+  close(conn)
+  
+  # fonction permettant de régler le problème de la délimitation entre colonnes (des fois un espace, des fois deux ou trois ...)
+  dummy.white.space = function(char.line){
+    temp = unlist(strsplit(char.line," ",fixed = TRUE))
+    test.temp = (temp =="")
+    result.vector = temp[!test.temp]
+    return(result.vector)
+  }
+  
+  start.line <- 6 # première ligne avec des données
+  end.line <- length(linn)-1 # le fichier finissant par un 'end', la dernière ligne de valeurs se situe juste à la ligne au-dessus.
+  nb.line <- end.line-(start.line-1) # nombre total de ligne
+  boundary.data <- data.frame(time = rep(0,nb.line), cum.seep.flux = rep(0,nb.line), cum.inj.flux = rep(0,nb.line), tot.seep.flux = rep(0,nb.line), tot.inj.flux = rep(0,nb.line), conc.seep = rep(0,nb.line), conc.inj = rep(0,nb.line)) # preallocation de la data.frame avec les résultats.
+  k <- 1 # compteur
+  # boucle sur les lignes avec des données
+  for (i in start.line:end.line){
+    current.line <- as.numeric(dummy.white.space(linn[i])) # transformation de la ligne en vecteur numérique
+    boundary.data$time[k] <- current.line[1] # time (h)
+    boundary.data$cum.seep.flux[k] <- -current.line[5] # ChemS2 [M L-1] cummulative solute flux across seepage face
+    boundary.data$cum.inj.flux[k] <- current.line[4] # ChemS3 [M L-1] cummulative solute flux injected
+    boundary.data$seep.flux[k] <- -current.line[3] # SMean2 [M T-1 L-1] total solute flux across seepage 
+    boundary.data$inj.flux[k] <- current.line[2] # SMean3 [M T-1 L-1] total solute flux injected
+    boundary.data$conc.seep[k] <- current.line[10] # cMean [M L-3] concentration seepage face
+    boundary.data$conc.inj[k] <- current.line[8] # cMean [M L-3] concentration injected 
+    k <- k+1 # incrementation du compteur
+  }
+  return(boundary.data)
+}
+
+recovered.ratio <- function(solute.data){
+  # injected mass
+  inj.mass <- solute.data$cum.inj.flux[dim(solute.data)[1]]
+  # recovered mass
+  recover.mass <- solute.data$cum.seep.flux[dim(solute.data)[1]]
+  recover.ratio <- recover.mass/inj.mass
+  return(recover.ratio)
+}
+
+plot.cum.break <- function(solute.data, t.inj, reference = "recovered"){
+  if (reference == "recovered"){
+    total.mass <- solute.data$cum.seep.flux[dim(solute.data)[1]]
+  } else if (reference == "injected") {
+    total.mass <- solute.data$cum.inj.flux[dim(solute.data)[1]]
+  } else {
+    stop("invalid argument value in reference")
+  }
+  plot(solute.data$time[solute.data$time > t.inj], solute.data$cum.seep.flux[solute.data$time > t.inj]/total.mass)
+}
+
+percentile.break <- function(solute.data, t.inj, level, reference = "recovered"){
+  if (reference == "recovered"){
+    total.mass <- solute.data$cum.seep.flux[dim(solute.data)[1]]
+  } else if (reference == "injected") {
+    total.mass <- solute.data$cum.inj.flux[dim(solute.data)[1]]
+  } else {
+    stop("invalid argument value in reference")
+  }
+  rel.cum.seep.flux <- solute.data$cum.seep.flux/total.mass
+  if (level > solute.data$cum.seep.flux[dim(solute.data)[1]]/total.mass){
+    rt <- NA
+  } else {
+    rt <- solute.data$time[which.min((rel.cum.seep.flux-level)^2)] - t.inj 
+  }
+  return(rt)
+}
+
+residence.time <- function(solute.data, t.inj){
+  inj.mass <- solute.data$cum.seep.flux[dim(solute.data)[1]]
+  rt.moy <- 1/inj.mass*trapz(solute.data$time,solute.data$time*solute.data$seep.flux) - t.inj
+  return(rt.moy)
+}
\ No newline at end of file
diff --git a/hydrus/processor_load.bat b/hydrus/processor_load.bat
new file mode 100644
index 0000000000000000000000000000000000000000..411b6d5a274697edd8167a2f5c23f851dea83080
--- /dev/null
+++ b/hydrus/processor_load.bat
@@ -0,0 +1 @@
+typeperf "\processeur(_Total)\%% temps processeur" -sc 5
diff --git a/hydrus/simulation_fcn.R b/hydrus/simulation_fcn.R
new file mode 100644
index 0000000000000000000000000000000000000000..fe28acc2bf3b628c3fed7c0157204b28ea2aeb48
--- /dev/null
+++ b/hydrus/simulation_fcn.R
@@ -0,0 +1,52 @@
+# This function process the output of the script .bat that estimates processor usage.
+output.bat <- function(output){
+  usage <- rep(0,5)
+  for (i in 1:5){
+    usage[i] <- as.numeric(strsplit(output[4+i],split='"',fixed=TRUE)[[1]][4])    
+  }
+  mean_usage <- mean(usage)
+  return(mean_usage)
+}
+
+# This function creates the file tree necessary to run simulation
+# It creates all simulation folders, copy files and modify the one containing the van Genuchten parameters
+simulation.setup <- function(nb.runs,generic.name,main.directory,hydrus.folder,new.data,hydrus.version="1D"){
+  for (i in 1:nb.runs){
+    # current simulation name
+    new.name <- paste(generic.name,as.character(i),sep = "")
+    # create a copy of the reference model in a folder named as the current simulation
+    new.hydrus.subfolder(new.name,main.directory,hydrus.folder)
+    # create the level_01 file
+    level_01.creator(new.name,main.directory)
+    # modification of the selector.in file
+    #new.data <- rbind(deposit,gravel)
+    modify.hydraulic.parameters(new.name,hydrus.folder,new.data[[i]],hydrus.version)
+  }
+}
+
+simulation.run <- function(nb.runs,processor_load_path,main_directory,generic.name,hydrus.calc="H1D_Calc.exe"){
+  start.time <- Sys.time()
+  sim.counter <- 0
+  
+  while (sim.counter < nb.runs){
+    # how busy the processor is?
+    temp <- system(processor_load_path,intern = TRUE)
+    processor.usage <- output.bat(temp)
+    if (processor.usage < 75){
+      # if processor usage is less than 75%, start a new run
+      sim.counter <- sim.counter + 1
+      # go to the simulation working directory
+      setwd(paste(main.directory,"/",generic.name,as.character(sim.counter),sep=""))
+      # run Hydrus
+      system(hydrus.calc, wait = FALSE)
+      # print the current simulation number
+      print(sim.counter)
+      # go back to main directory
+      setwd(main.directory)
+    }
+    # wait a minute before restarting the loop
+    Sys.sleep(60)
+  }
+  end.time <- Sys.time()
+  return(end.time-start.time)
+}
\ No newline at end of file