vignettes/boundary_conditions.Rmd
boundary_conditions.Rmd##
! {RMODFLOW} is still in its experimental lifecycle stage.
##
! Use at your own risk, and submit issues here:
##
! <https://github.com/rogiersbart/RMODFLOW/issues>Boundary conditions are a inherent part of partial differential equations. In terms of groundwater flow modelling, they typically represent features that add fluxes to the system or specify constant head values. These features may vary in space and time. MODFLOW is quite flexible in the way it lets users define boundary conditions. There are multiple ways to ‘group’ boundary conditions (e.g. Dirichlet, Neumann, Cauchy types). For the purpose of the RMODFLOW API however, all MODFLOW boundary conditions can be placed in one of two groups:
This is similar to the way RMODFLOW handles input data (see also the vignette in input data structures), i.e. rmf_lists and rmf_arrays.
Some of the basic boundary condition packages in MODFLOW-2005 are (for an extensive overview of RMODFLOW supported packages, check out https://github.com/rogiersbart/RMODFLOW/wiki/Overview-of-supported-MODFLOW-packages):
The RCH & EVT packages represent spatially continuous data; the others discrete features.
As with other packages, RMODFLOW can read in boundary condition packages with the appropriate rmf_read_* function:
dis <- rmf_example_file('water-supply-problem.dis') %>%
rmf_read_dis()
riv <- rmf_example_file('water-supply-problem.riv') %>%
rmf_read_riv(dis = dis)The discrete boundary condition packages are lists with following elements:
dimensions holding information on the amount of features & parametersoption with possible optionsaux with the names of possible AUX variablesdata a rmf_list with the data of the boundary condition featureskper a data.frame with the stress-period information of the features
str(riv)
#> List of 10
#> $ np : num 0
#> $ mxl : num 0
#> $ instances: NULL
#> $ mxact : num 67
#> $ itmp : Named num 67
#> ..- attr(*, "names")= chr "list_1"
#> $ irivcb : num 9
#> $ option : Named logi FALSE
#> ..- attr(*, "names")= chr "noprint"
#> $ aux : chr "IFACE"
#> $ data :Classes 'rmf_list' and 'data.frame': 67 obs. of 9 variables:
#> ..$ k : int [1:67] 1 1 1 1 1 1 1 1 1 1 ...
#> ..$ i : int [1:67] 18 19 19 19 19 20 20 20 21 21 ...
#> ..$ j : int [1:67] 1 1 2 3 4 4 5 6 6 7 ...
#> ..$ stage : num [1:67] 0 0 0 0 0 0 0 0 0 0 ...
#> ..$ conductance: num [1:67] 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 ...
#> ..$ rbot : num [1:67] -10 -10 -10 -10 -10 -10 -10 -10 -10 -10 ...
#> ..$ IFACE : num [1:67] 0 0 0 0 0 0 0 0 0 0 ...
#> ..$ parameter : logi [1:67] FALSE FALSE FALSE FALSE FALSE FALSE ...
#> ..$ name : chr [1:67] "list_1" "list_1" "list_1" "list_1" ...
#> $ kper :'data.frame': 1 obs. of 2 variables:
#> ..$ kper : int 1
#> ..$ list_1: logi TRUE
#> - attr(*, "class")= chr [1:2] "riv" "rmf_package"
#> - attr(*, "comment")= chr " RIV: River package file created on 2017-05-30 by ModelMuse version 3.9.0.0."Writing can be done with the appropriate rmf_write_* function:
rmf_write_riv(riv, dis = dis, file = 'input.riv')Creating discrete boundary condition objects can be done by supplying rmf_list objects with the appropriate variables and kper argument to the rmf_create_* function:
riv_east <- data.frame(i = 1:dis$nrow, j = dis$ncol, k = 1, stage = 0, conductance = 0.02, rbot = -10) %>%
rmf_create_list(kper = 1)
riv_west <- data.frame(i = 1:dis$nrow, j = 1, k = 1, stage = 0, conductance = 0.02, rbot = -10) %>%
rmf_create_list(kper = 1)
riv_new <- rmf_create_riv(riv_east, riv_west, dis = dis)
str(riv_new)
#> List of 9
#> $ np : num 0
#> $ mxl : num 0
#> $ instances: NULL
#> $ mxact : num 90
#> $ itmp : Named num [1:2] 45 45
#> ..- attr(*, "names")= chr [1:2] "list_1" "list_2"
#> $ irivcb : num 0
#> $ option : Named logi FALSE
#> ..- attr(*, "names")= chr "noprint"
#> $ data :Classes 'rmf_list' and 'data.frame': 90 obs. of 8 variables:
#> ..$ k : num [1:90] 1 1 1 1 1 1 1 1 1 1 ...
#> ..$ i : int [1:90] 1 2 3 4 5 6 7 8 9 10 ...
#> ..$ j : num [1:90] 56 56 56 56 56 56 56 56 56 56 ...
#> ..$ stage : num [1:90] 0 0 0 0 0 0 0 0 0 0 ...
#> ..$ conductance: num [1:90] 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 ...
#> ..$ rbot : num [1:90] -10 -10 -10 -10 -10 -10 -10 -10 -10 -10 ...
#> ..$ parameter : logi [1:90] FALSE FALSE FALSE FALSE FALSE FALSE ...
#> ..$ name : chr [1:90] "list_1" "list_1" "list_1" "list_1" ...
#> $ kper :'data.frame': 1 obs. of 3 variables:
#> ..$ kper : int 1
#> ..$ list_1: logi TRUE
#> ..$ list_2: logi TRUE
#> - attr(*, "class")= chr [1:2] "riv" "rmf_package"As input, the rmf_create_* functions can either take one or more rmf_list objects (possibly of class rmf_parameter) or a single list with one or more rmf_list objects (possibly of class rmf_parameter). If no kper argument is set on any of the objects, RMODFLOW will throw an error:
riv_east <- data.frame(i = 1:dis$nrow, j = dis$ncol, k = 1, stage = 0, conductance = 0.02, rbot = -10) %>%
rmf_create_list()
riv_new <- rmf_create_riv(riv_east, dis = dis)
#> Error: Please make sure all rmf_list and rmf_parameter objects have a kper attributeSimilar to their discrete kin, continuous boundary condition packages can be read in through the appropriate rmf_read_* function and written with their rmf_write_* function. When reading a continuous boundary condition package that has been defined using MODFLOW parameters (see also the vignette on RMODFLOW parameters), a mlt and zon argument should often be specified.
# Currently, no example models with RCH or EVT packagesCreating continuous boundary condition packages can be done by supplying rmf_2d_array objects (possibly of class rmf_parameter) or a single list with rmf_2d_array objects (possibly of class rmf_parameter) to the rmf_create_* function. Similar to the discrete boundary conditions, a kper argument must be set on all objects:
dis <- rmf_create_dis(nper = 2)
rch_summer <- rmf_create_array(0.0002, dim = c(dis$nrow, dis$ncol), kper = 1)
rch_winter <- rmf_create_array(0.0003, dim = c(dis$nrow, dis$ncol), kper = 2)
rch <- rmf_create_rch(list(rch_summer, rch_winter), dis = dis)The continuous boundary condition packages are lists with following elements:
dimensions holding information on the amount of features & parametersrecharge which is a list with rmf_2d_arrays holding the recharge valuesirch which are also lists with rmf_2d_arrays holding those valueskper a data.frame with the stress-period information of the features
str(rch)
#> List of 7
#> $ np : num 0
#> $ instances: NULL
#> $ nrchop : num 3
#> $ irchcb : num 0
#> $ recharge :List of 2
#> ..$ array_1: 'rmf_2d_array' num [1:10, 1:10] 2e-04 2e-04 2e-04 2e-04 2e-04 2e-04 2e-04 2e-04 2e-04 2e-04 ...
#> .. ..- attr(*, "kper")= num 1
#> .. ..- attr(*, "dimlabels")= chr [1:2] "i" "j"
#> ..$ array_2: 'rmf_2d_array' num [1:10, 1:10] 3e-04 3e-04 3e-04 3e-04 3e-04 3e-04 3e-04 3e-04 3e-04 3e-04 ...
#> .. ..- attr(*, "kper")= num 2
#> .. ..- attr(*, "dimlabels")= chr [1:2] "i" "j"
#> $ kper :'data.frame': 2 obs. of 3 variables:
#> ..$ kper : int [1:2] 1 2
#> ..$ array_1: logi [1:2] TRUE FALSE
#> ..$ array_2: logi [1:2] FALSE TRUE
#> $ irchpf : num -1
#> - attr(*, "class")= chr [1:2] "rch" "rmf_package"Notice that for different stress-periods, the variables (e.g. recharge) are defined as separate rmf_2d_arrays with kper attributes rather than a single 3D array where the 3th dimension represents stress-periods.