Groundwater flow simulation and particle tracking to forecast microbial concentration in a drinking water extraction well.
- "gdal_polygonize.py" must be installed (-->
install python-gdal
) - MODFLOW and MODPATH must be installed (
mfusg
for MODFLOW usgs,mp6
for MODPATH). If you need binaries compiled for linux (ubuntu) please contact me - R must be installed
Two possiblities:
-
Open R, copy-paste
run.R
or sourcerun.R
(but don't forget to adapt line 15 inrun.R
theDIR
variable to your directory structure -
Open terminal, change directory to "gwModBac" (
cd path/gwModBac
), enterRscript --vanilla run.R sim_name 15 50 10 4 NULL 100 200 20 OK
run.R
--> path of the filerun.R
sim_name
name of the simulation (an output filesim_name.txt
will be created in the directorygwModBac
)15
--> nx: number of cells along the x-direction = number of columns50
--> ny: number of cells along the y-direction = number of rows10
--> nz: number of cells along the z-direction = number of layers4
--> number of runs (but one output file with all the output together)NULL
--> no seed are used. To set a seed for the random number simulator, replaceNULL
by any numeric value (e.g.,10
). With a seed, you can reproduce the simulation exactly.100
--> nx_ref: number of columns of the reference simulation (nx_ref >= nx), used to scale down the Gaussian Random Field (only used if a seed number is set)200
--> ny_ref: number of rows of the reference simulation (ny_ref >= ny), used to scale down the Gaussian Random Field (only used if a seed number is set)20
--> nz_ref: number of layers of the reference simulation (nz_ref >= nz), used to scale down the Gaussian Random Field (only used if a seed number is set)OK
[optional] --> can be anything; if the 7th argument is present, an overview plot is created for each single simulation (in the directorysimulations
): same plot as in section Groundwater flow simulation and particle tracking. If absent, no plot will be created (faster).
Example with fixed seed:
```
Rscript --vanilla run.R sim01 10 20 10 1 102343643 100 200 10 OK
```
Example without fixed seed and without plot:
```
Rscript --vanilla run.R sim01 10 20 10 1 NULL 100 200 10
```
Only one argument: the config file *.yaml
. All the parameters as well as their prior distributions are explained in the test file sim_parameters.yaml
.
To simulate the 3D hydraulic conductivity field (random Gaussian field), you have two options:
-
either assign the values the 3D random Gaussian field parameter (
para: TRUE
) -
or give directly the values of the 3D random Gaussian field parameter (
para: FALSE
and the code read the file given byfile: 'HK.txt'
): one big vector column of length (nx x ny x nz), corresponding to an array of dimension ny x nx x nz ( = nrow x ncol x nlayers, notice that here nrow = ny). I haven't tested what is the order of this big vector.Rscript --vanilla run_para.R sim_parameters.yaml
Output:
- a text file with name
projectName/simName.txt
(seesim_parameters.yaml
) with a vector column corresponding to the 10-days ahead forecast - if
plot: TRUE
a plotprojectName/simName.png
Minimal possible grid size:
- nx >= 10
- ny >= 10
- nz >= 2
- Read the hyperparameters/parameters from
para.R
- Read the data (river stage time-series, groundwater head time-series) from
data/
- Create the simulation grid (according to parameters)
- Define
- the river properties (location, stage, riverbed),
- the location of the specified head boundary,
- the particles to estimate the microbila concentration in the drinking water extraction well.
- Run Monte Carlo simulations (unconditional)
- simulate the hydraulic properties: hydraulic conductivity (Gaussian random field), porosity, etc.
- simulate the specified head boundary conditions (Gaussian process)
for the 10-days ahead forecast:
- simulate the precipitation for the next 10 days
- using a convolution model between precipitation and river stage, simulate the river stage for the next 10 days
- using a convolution model between river stage and groundwater heads at the observation wells, simulate the groundwater heads at the observation wells for the next 10 days
- Simulate the head boundary conditions conditional on the groundwater heads the observation wells
- Interpolate from the head boundary conditions the initial heads
- run MODFLOW
- run MODPATH to simulate pathway of the microbes that reach the well drinking water extraction well
- Simulate the microbial concentration in the river from the river stage
- Predict the microbial concentration in the well according to an exponential decay model that simulates the mechanical filtration of microbes in the aquifer
- contours and colour bar = hydraulic heads of the first layer
- coloured lines = paths of the particles entering into the well
- green dots = start positions of the particles coming from the river
- red dots = start positions of the particles that do not come from the river
- light blue multipolygon (corner at (0, 0), (0, 150), (5, 150), (5, 0) = the river
Some simulation results (with seed = 102343643):
- Target value as a function of nx, ny, nz:
- Target value as a function of nx, ny for nz = 10 (z-axis = target value!)
For these results, I recommand to keep nz = 10 and to only vary nx and ny (more stable, better convergence).