## 
! {RMODFLOW} is still in its experimental lifecycle stage.
## 
! Use at your own risk, and submit issues here:
## 
! <https://github.com/rogiersbart/RMODFLOW/issues>

Since MODFLOW-2000, you can specify so called parameters in the MODFLOW input files of flow and boundary condition packages. A parameter multiplies all the variable values of the defined cells with the parameter value. For example in the river package, the only allowed parameter type is RIV, which represents the riverbed hydraulic conductivity. All the riverbed hydraulic conductivities in the cells defined under the parameter will therefore be multiplied with the parameter value. You can specify parameters for the horizontal hydraulic conductivity in layers, for the recharge values, the riverbed conductance etc. MODFLOW parameters differ between discrete boundary conditions and continuous data for boundary conditions and flow packages. These parameters all have the following:

  • A parameter name
  • A predefined parameter type specific to the package (e.g. HK, VK, RCH, …)
  • A parameter value
  • A number specifying the amount of clusters/cells defined by the parameter

Additionally, a parameter can be time-varying. This means that there can be multiple “instances” of a single parameter. Each instance defines a different set of clusters/cells for the parameter. The name, type and value remain the same.

The idea behind MODFLOW parameters is that they are easier to handle during calibration and sensitivity analysis. Instead of changing all individual values of each cell, you just change the parameter value. This is made even easier by the PVAL file, which lists all the parameters active in the current simulation with their names and values. This file can be used by external optimization softwares such as UCODE.

If you do not use a scripting language to prepare or optimize your MODFLOW model, parameters can be very useful. But if you use a scripting language you could just as easily define a parameter inside your script and use those in your scripted optimization.

I strongly advise against using MODFLOW parameters in your models, unless you really need to. You can script them yourself and keep your input files self-contained.

MODFLOW parameters are difficult to handle because sometimes, the input specifying a single parameter is spread around multiple files. With a single function, it is difficult to read in or access a single MODFLOW parameter. Your new best friend RMODFLOW is equipped to deal with this however.

The RMODFLOW parameter API aims at :

  1. Creating objects of class rmf_parameter that are self-contained, i.e. all the data required by the parameter is contained in a single R object. This is done by setting attributes on the object such as parnam, parval, mlt etc. This way, if a user wants to access a MODFLOW parameter, they don’t have to combine or calculate from different objects.
  2. Creating these rmf_parameter objects with RMODFLOW using S3 functions which should be relatively straightforward (or maybe not, as we’ll see below).

Reading

Reading of MODFLOW parameters is done automatically by the RMODFLOW rmf_read_* of those packages that may contain parameters.

For discrete boundary conditions, there is no problem. A parameter is defined with data much like it is without a parameter. The data object also contains a logical parameter column indicating if the cell is part of a parameter. A new parameter_values element is also added.

dis <- rmf_example_file("example-model.dis") %>%
  rmf_read_dis()

# chd <- rmf_example_file("example-model.chd") %>%
#   rmf_read_chd(dis = dis) # FIX failing for me at the moment
# str(chd)

The problem arises when reading/writing continuous data defined by parameters, i.e. flow package data and continuous boundary condition data. Parameters of this kind are defined for clusters, not cells. There can be multiple clusters per parameter. The information for a cluster is different for flow packages and continuous boundary conditions.

For flow packages, a parameter cluster contains:

  • Integers specifying the layer(s) to which the parameter applies
  • Name of multiplier array (“NONE” if no multiplier array is active)
  • Name of zone array (“ALL” if no zone array is active)
  • IZ values which are the integer indices of the zone array, if present, for which the parameter applies

For a continuous boundary condition, a parameter cluster contains:

  • Name of multiplier array (“NONE” if no multiplier array is active)
  • Name of zone array (“ALL” if no zone array is active)
  • IZ values which are the integer indices of the zone array, if present, for which the parameter applies
  • Additionally, this type of parameter can be time-varying.s Note: time-varying parameters are implemented but have not been tested.

The problem here is that if we want to read in a package defined by such parameters, we need to supply a mlt and/or zon object to the read function (unless mltarr = “NONE” and zonarr = “ALL”). The idea is that the read-in object contains ALL the information used by MODFLOW. So not only the multiplier/zone array names, but also their values. That way, the user can access the parameter directly and in its entirety without having to calculate it from different arrays.

mlt <- rmf_example_file("example-model.mlt") %>%
  rmf_read_mlt(dis = dis)

zon <- rmf_example_file("example-model.zon") %>%
  rmf_read_zon(dis = dis)

lpf <- rmf_example_file("example-model.lpf") %>%
  rmf_read_lpf(dis = dis, mlt = mlt, zon = zon)
str(lpf)
#> List of 19
#>  $ ilpfcb            : num 9
#>  $ hdry              : num -2e+20
#>  $ nplpf             : num 2
#>  $ storagecoefficient: logi FALSE
#>  $ constantcv        : logi FALSE
#>  $ thickstrt         : logi FALSE
#>  $ nocvcorrection    : logi FALSE
#>  $ novfc             : logi FALSE
#>  $ noparcheck        : logi FALSE
#>  $ laytyp            : num [1:9] 0 0 0 0 0 0 0 0 0
#>  $ layavg            : num [1:9] 0 0 0 0 0 0 0 0 0
#>  $ chani             : num [1:9] -1 -1 -1 -1 -1 -1 -1 -1 -1
#>  $ layvka            : num [1:9] 0 0 0 0 0 0 0 0 0
#>  $ laywet            : num [1:9] 0 0 0 0 0 0 0 0 0
#>  $ parameters        :List of 2
#>   ..$ HK_Par1: 'rmf_parameter' num [1:10, 1:10, 1:6] 0 0 0 0 0 0 0 0 0 0 ...
#>   .. ..- attr(*, "dimlabels")= chr [1:3] "i" "j" "k"
#>   .. ..- attr(*, "parnam")= chr "HK_Par1"
#>   .. ..- attr(*, "parval")= num 150
#>   .. ..- attr(*, "layer")= num [1:6] 1 2 3 4 5 6
#>   .. ..- attr(*, "partyp")= chr "HK"
#>   .. ..- attr(*, "mlt")= chr [1:6] "NONE" "NONE" "NONE" "NONE" ...
#>   .. ..- attr(*, "zon")= chr [1:6] "HK_Par1_Z1" "HK_Par1_Z2" "HK_Par1_Z3" "HK_Par1_Z4" ...
#>   .. ..- attr(*, "iz")=List of 6
#>   .. .. ..$ : num 1
#>   .. .. ..$ : num 1
#>   .. .. ..$ : num 1
#>   .. .. ..$ : num 1
#>   .. .. ..$ : num 1
#>   .. .. ..$ : num 1
#>   ..$ HK_Par2: 'rmf_parameter' num [1:10, 1:10, 1:9] 100 100 100 100 90 ...
#>   .. ..- attr(*, "dimlabels")= chr [1:3] "i" "j" "k"
#>   .. ..- attr(*, "parnam")= chr "HK_Par2"
#>   .. ..- attr(*, "parval")= num 200
#>   .. ..- attr(*, "layer")= num [1:9] 1 2 3 4 5 6 7 8 9
#>   .. ..- attr(*, "partyp")= chr "HK"
#>   .. ..- attr(*, "mlt")= chr [1:9] "HK_Par2_M1" "HK_Par2_M2" "HK_Par2_M3" "HK_Par2_M4" ...
#>   .. ..- attr(*, "zon")= chr [1:9] "HK_Par2_Z1" "HK_Par2_Z2" "HK_Par2_Z3" "HK_Par2_Z4" ...
#>   .. ..- attr(*, "iz")=List of 6
#>   .. .. ..$ : num 1
#>   .. .. ..$ : num 1
#>   .. .. ..$ : num 1
#>   .. .. ..$ : num 1
#>   .. .. ..$ : num 1
#>   .. .. ..$ : num 1
#>  $ parameter_values  : Named num [1:2] 150 200
#>   ..- attr(*, "names")= chr [1:2] "HK_Par1" "HK_Par2"
#>  $ hk                : 'rmf_3d_array' num [1:10, 1:10, 1:9] 100 100 100 100 90 ...
#>   ..- attr(*, "dimlabels")= chr [1:3] "i" "j" "k"
#>  $ hani              : 'rmf_3d_array' num [1:10, 1:10, 1:9] 1 1 1 1 1.11 ...
#>   ..- attr(*, "dimlabels")= chr [1:3] "i" "j" "k"
#>  $ vka               : 'rmf_3d_array' num [1:10, 1:10, 1:9] 1 1 1 1 1 1 1 1 1 1 ...
#>   ..- attr(*, "dimlabels")= chr [1:3] "i" "j" "k"
#>  - attr(*, "comment")= chr " LPF: Layer Property Flow package file created on 2017-05-30 by ModelMuse version 3.9.0.0."
#>  - attr(*, "class")= chr [1:2] "lpf" "rmf_package"

If a mlt or zon object is needed but not supplied, the rmf_read_* function will throw an error.

The returned flow package object will have a parameter element which contains the flow parameters and a parameter_values argument containing the corresponding values. Additionally, the variable (let’s call them variables here; there actually hydraulic parameters since they remain constant during the simulation) arrays such as HK have been recalculated taking the defined parameters into account. Thus, in the example above, where parameters were defined for HK, the lpf$hk array contains the values used by MODFLOW. You don’t have to calculate those yourself.

For continuous boundary condition packages such as recharge, the returned object also contains a parameter_values element but the parameter arrays are listed under the corresponding variable, e.g. rch$recharge. When plotting, the S3 method automatically calculates and plots all the active arrays for the specified stress-period.

# No example with continuous boundary package defined by parameters

If you do not want to deal with the reading of additional input arguments, just use the top-level rmf_read() function as discussed in the vignette on top-level functions.

You can calculate a resulting rmf_array from specified multiplier & zone arrays, parameter values and or layers indices using rmf_calculate_array():

rmf_calculate_array(dis = dis, 
                    layer = 1:2,
                    mltarr = mlt$rmlt[c(1, 9)],
                    zonarr = NULL,
                    parval = 1.5)
#> RMODFLOW 3d array with 10 rows, 10 columns and 2 layers, representing the i, j & k dimensions. 
#> Not representing stress period data 
#> , , 1
#> 
#>            [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
#>  [1,] 0.7500000 0.7500000 0.7500000 0.7500000 0.7686140 0.8474689 0.9263238
#>  [2,] 0.7500000 0.7500000 0.7500000 0.7500000 0.7893653 0.8682202 0.9470751
#>  [3,] 0.7500000 0.7500000 0.7500000 0.7500000 0.8101166 0.8889715 0.9678264
#>  [4,] 0.7500000 0.7432230 0.7222131 0.7012033 0.7815875 0.8979043 0.9885777
#>  [5,] 0.6749410 0.6539311 0.6329213 0.6119115 0.6599835 0.7763003 0.8942276
#>  [6,] 0.5856492 0.5646393 0.5436295 0.5226197 0.5383795 0.6597886 0.7900325
#>  [7,] 0.4963574 0.4753475 0.4543377 0.4500000 0.4500000 0.5555935 0.6858374
#>  [8,] 0.4500000 0.4500000 0.4500000 0.4500000 0.4500000 0.4513984 0.5816423
#>  [9,] 0.4500000 0.4500000 0.4500000 0.4500000 0.4500000 0.4500000 0.4774472
#> [10,] 0.4500000 0.4500000 0.4500000 0.4500000 0.4500000 0.4500000 0.4500000
#>            [,8]      [,9]     [,10]
#>  [1,] 1.0051788 1.0500000 1.0500000
#>  [2,] 1.0259301 1.0500000 1.0500000
#>  [3,] 1.0466813 1.0500000 1.0500000
#>  [4,] 1.0500000 1.0500000 1.0500000
#>  [5,] 1.0244715 1.0500000 1.0500000
#>  [6,] 0.9202764 1.0500000 1.0500000
#>  [7,] 0.8160813 0.9463252 1.0500000
#>  [8,] 0.7118862 0.8421301 0.9723740
#>  [9,] 0.6076911 0.7379350 0.8681789
#> [10,] 0.5034959 0.6337398 0.7639837
#> 
#> , , 2
#> 
#>            [,1]      [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
#>  [1,] 0.7500000 0.7500000 0.7500000 0.7500000 0.7686140 0.8474689 0.9263238
#>  [2,] 0.7500000 0.7500000 0.7500000 0.7500000 0.7893653 0.8682202 0.9470751
#>  [3,] 0.7500000 0.7500000 0.7500000 0.7500000 0.8101166 0.8889715 0.9678264
#>  [4,] 0.7500000 0.7432230 0.7222131 0.7012033 0.7815875 0.8979043 0.9885777
#>  [5,] 0.6749410 0.6539311 0.6329213 0.6119115 0.6599835 0.7763003 0.8942276
#>  [6,] 0.5856492 0.5646393 0.5436295 0.5226197 0.5383795 0.6597886 0.7900325
#>  [7,] 0.4963574 0.4753475 0.4543377 0.4500000 0.4500000 0.5555935 0.6858374
#>  [8,] 0.4500000 0.4500000 0.4500000 0.4500000 0.4500000 0.4513984 0.5816423
#>  [9,] 0.4500000 0.4500000 0.4500000 0.4500000 0.4500000 0.4500000 0.4774472
#> [10,] 0.4500000 0.4500000 0.4500000 0.4500000 0.4500000 0.4500000 0.4500000
#>            [,8]      [,9]     [,10]
#>  [1,] 1.0051788 1.0500000 1.0500000
#>  [2,] 1.0259301 1.0500000 1.0500000
#>  [3,] 1.0466813 1.0500000 1.0500000
#>  [4,] 1.0500000 1.0500000 1.0500000
#>  [5,] 1.0244715 1.0500000 1.0500000
#>  [6,] 0.9202764 1.0500000 1.0500000
#>  [7,] 0.8160813 0.9463252 1.0500000
#>  [8,] 0.7118862 0.8421301 0.9723740
#>  [9,] 0.6076911 0.7379350 0.8681789
#> [10,] 0.5034959 0.6337398 0.7639837

As stated above, a PVAL file contains all the active MODFLOW parameters. This can be read using the rmf_read_pval() function:

pvl <- rmf_example_file("example-model.pval") %>%
  rmf_read_pval()
str(pvl)
#> List of 2
#>  $ np  : num 4
#>  $ data: tibble [4 x 2] (S3: tbl_df/tbl/data.frame)
#>   ..$ parnam: chr [1:4] "HK_Par1" "HK_Par2" "CHD_Par1" "CHD_Par2"
#>   ..$ parval: num [1:4] 150 200 30 30
#>  - attr(*, "comment")= chr " PVAL file created on 2017-05-30 by ModelMuse version 3.9.0.0."
#>  - attr(*, "class")= chr [1:2] "pval" "rmf_package"

Writing

Once the object is read in or created, writing is no problem. Just use the appropriate rmf_write_* function and this will handle everything. No need to supply additional function arguments besides the dis object and the filepath.

Creating

RMODFLOW provides S3 methods for creating MODFLOW parameters through the rmf_create_parameter() functions.

Creating a parameter for a discrete boundary condition package is easy: just call rmf_create_parameter.rmf_list() on a rmf_list object:

lst <- data.frame(i = 1,
                  j = 1:dis$ncol,
                  k = 1,
                  q = 500) %>%
  rmf_create_list(kper = 1)

lst_p <- rmf_create_parameter(lst, parnam = "WELLS", parval = 1.5)
str(lst_p)
#> Classes 'rmf_parameter', 'rmf_list' and 'data.frame':    10 obs. of  4 variables:
#>  $ i: num  1 1 1 1 1 1 1 1 1 1
#>  $ j: int  1 2 3 4 5 6 7 8 9 10
#>  $ k: num  1 1 1 1 1 1 1 1 1 1
#>  $ q: num  500 500 500 500 500 500 500 500 500 500
#>  - attr(*, "kper")= num 1
#>  - attr(*, "parnam")= chr "WELLS"
#>  - attr(*, "parval")= num 1.5

Creating a parameter for continuous boundary conditions is a little trickier. There are two ways of doing this. First, you can call rmf_create_parameter.rmf_2d_array on a rmf_2d_array and specify the kper argument which by default looks at the kper argument of the array. This is intended to be used on arrays that contain the same values. If values differ in an array however, the mltnam attribute is named (which can be set with the mltnam argument). As a user, you should then supply this array to a MLT object.

ar <- rmf_create_array(12e-1, dim = c(dis$nrow, dis$ncol)) %>%
  rmf_create_parameter(parnam = "RCH_1", parval = 1.5, kper = 1) 
ar
#> RMODFLOW 2d array with 10 rows and 10 columns, representing the i & j dimensions. 
#> Active during stress period: 1 
#>       [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#>  [1,]  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2   1.2
#>  [2,]  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2   1.2
#>  [3,]  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2   1.2
#>  [4,]  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2   1.2
#>  [5,]  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2   1.2
#>  [6,]  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2   1.2
#>  [7,]  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2   1.2
#>  [8,]  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2   1.2
#>  [9,]  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2   1.2
#> [10,]  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2  1.2   1.2
#> attr(,"parnam")
#> [1] "RCH_1"
#> attr(,"parval")
#> [1] 1.5
#> attr(,"mlt")
#> [1] "NONE"
#> attr(,"zon")
#> [1] "ALL"

ar2 <- rnorm(dis$nrow*dis$ncol, mean = 12e-1, sd = 0.5) %>%
  rmf_create_array(dim = c(dis$nrow, dis$ncol)) %>%
  rmf_create_parameter(parnam = "RCH_2", parval = 1.5, kper = 1) 
ar2
#> RMODFLOW 2d array with 10 rows and 10 columns, representing the i & j dimensions. 
#> Active during stress period: 1 
#>            [,1]       [,2]      [,3]      [,4]      [,5]      [,6]      [,7]
#>  [1,] 1.3987215  1.4007822 0.9481330 1.1861319 0.5285347 1.0197793 1.3316558
#>  [2,] 1.1014303  0.3172325 1.2282764 0.7300679 1.6681879 1.1518832 0.7751159
#>  [3,] 1.6531710  0.3876060 1.2060965 1.0574406 1.2965436 0.9930334 1.4943110
#>  [4,] 0.7814193  1.5404219 1.8328985 0.9367031 0.8126319 0.9844248 1.4360198
#>  [5,] 1.2394528  1.6142921 1.5891760 1.3831091 2.2647125 1.3270520 0.7856721
#>  [6,] 1.7520531  1.0204952 1.2735220 1.6790536 0.5510337 1.6078179 1.4659023
#>  [7,] 1.3560670  0.7346869 1.0309370 1.4310345 1.3094746 2.0020454 0.6993674
#>  [8,] 1.3241762  0.7967068 1.2289476 1.2881479 1.5803410 1.1627943 1.6313164
#>  [9,] 1.0465408  0.6020116 0.8228026 1.0803232 1.1579018 0.5282616 1.0703940
#> [10,] 0.9539094 -0.1537345 1.0997940 0.9129385 0.6730613 1.5048293 1.9936104
#>            [,8]     [,9]     [,10]
#>  [1,] 1.2530623 1.352850 1.3610880
#>  [2,] 0.8609448 1.228679 1.7992211
#>  [3,] 2.0252072 1.277916 1.6909336
#>  [4,] 0.2547936 0.665999 0.3178983
#>  [5,] 1.1027948 0.809356 1.7271896
#>  [6,] 1.2095580 1.347287 1.1084517
#>  [7,] 0.5601449 1.371613 0.9257086
#>  [8,] 1.6866873 1.454696 0.1696303
#>  [9,] 0.5172566 1.669242 0.3597518
#> [10,] 1.5336927 1.123371 0.9521509
#> attr(,"parnam")
#> [1] "RCH_2"
#> attr(,"parval")
#> [1] 1.5
#> attr(,"mlt")
#> [1] "RCH_2"
#> attr(,"zon")
#> [1] "ALL"

Secondly, you can use rmf_create_parameter.default() and supply the multiplier and/or zone arrays directly to create an array from:

mlt <- rmf_create_mlt(mltnam = "RCH", rmlt = rmf_create_array(0.0002, dim = c(dis$nrow, dis$ncol)))
zon <- rmf_create_zon(zonnam = "RCH", izon = rmf_create_array(1:4, dim = c(dis$nrow, dis$ncol)))
p <- rmf_create_parameter(dis = dis, kper = 1, parnam = "RCH_1", parval = 1.5,
                          mltnam = "RCH", zonnam = "RCH", iz = list(c(1,2)),
                          mlt = mlt, zon = zon)
p
#> RMODFLOW 2d array with 10 rows and 10 columns, representing the i & j dimensions. 
#> Active during stress period: 1 
#>        [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
#>  [1,] 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00
#>  [2,] 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00
#>  [3,] 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04
#>  [4,] 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04
#>  [5,] 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00
#>  [6,] 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00
#>  [7,] 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04
#>  [8,] 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04
#>  [9,] 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00
#> [10,] 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00 3e-04 0e+00
#> attr(,"parnam")
#> [1] "RCH_1"
#> attr(,"parval")
#> [1] 1.5
#> attr(,"mlt")
#> [1] "RCH"
#> attr(,"zon")
#> [1] "RCH"
#> attr(,"iz")
#> attr(,"iz")[[1]]
#> [1] 1 2

For flow package data, it works similar. In addition to the above functions, you also have rmf_create_parameter.rmf_3d_array(). Also, for flow package data, you should not supply kper but rather a layer argument (defaults to all layers in the array) as well as a partyp argument. Just like rmf_create_parameter.rmf_2d_array, this is intended to be used on arrays the contain the same values so that Mltarr = “NONE” and Zonarr = “ALL”. If different values are present, you should remember to add the array to the mlt object.

hk <- rmf_create_array(10, dim = c(dis$nrow, dis$ncol, dis$nlay)) %>%
  rmf_create_parameter(parnam = "HK_1", parval = 0.5, layer = 1:dis$nlay,
                       partyp = "HK")
str(hk)
#>  'rmf_parameter' num [1:10, 1:10, 1:9] 10 10 10 10 10 10 10 10 10 10 ...
#>  - attr(*, "dimlabels")= chr [1:3] "i" "j" "k"
#>  - attr(*, "parnam")= chr "HK_1"
#>  - attr(*, "parval")= num 0.5
#>  - attr(*, "layer")= int [1:9] 1 2 3 4 5 6 7 8 9
#>  - attr(*, "partyp")= chr "HK"
#>  - attr(*, "mlt")= chr [1:9] "NONE" "NONE" "NONE" "NONE" ...
#>  - attr(*, "zon")= chr [1:9] "ALL" "ALL" "ALL" "ALL" ...

It is therefore often safer to use the rmf_create_parameter.default() function:

mlt <- rmf_create_mlt(mltnam = "HK", rmlt = rmf_create_array(rnorm(dis$nrow*dis$ncol, 8, 8), dim = c(dis$nrow, dis$ncol)))
zon <- rmf_create_zon(zonnam = "HK", izon = rmf_create_array(1, dim = c(dis$nrow, dis$ncol)))
p <- rmf_create_parameter(dis = dis, layer = 1:dis$nlay, parnam = "HK_1", parval = 1.5,
                          mltnam = "HK", zonnam = "HK", iz = list(1),
                          partyp = "HK", mlt = mlt, zon = zon)
str(p)
#>  'rmf_parameter' num [1:10, 1:10, 1] -2.73 9.35 7.91 6.68 12.53 ...
#>  - attr(*, "dimlabels")= chr [1:3] "i" "j" "k"
#>  - attr(*, "parnam")= chr "HK_1"
#>  - attr(*, "parval")= num 1.5
#>  - attr(*, "layer")= int [1:9] 1 2 3 4 5 6 7 8 9
#>  - attr(*, "partyp")= chr "HK"
#>  - attr(*, "mlt")= chr "HK"
#>  - attr(*, "zon")= chr "HK"
#>  - attr(*, "iz")=List of 1
#>   ..$ : num 1

To create a pvl object, simply supply vectors with the names and parameter values:

rmf_create_pval(parnam = c("HK_1", "HK_2"),
               parval = c(12, 0.01))
#> RMODFLOW Parameter Value File object with: 
#> 2 parameter values: 
#> 
#>   parnam parval
#> 1   HK_1  12.00
#> 2   HK_2   0.01

Of course, all of this can be avoided by not using MODFLOW parameters and just script your own:

hk_fact <- 1.5
hk <- rmf_create_array(10, dim = c(dis$nrow, dis$ncol, dis$nlay)) * hk_fact
lpf <- rmf_create_lpf(dis = dis, hk = hk)