pypest / pyemu

python modules for model-independent uncertainty analyses, data-worth analyses, and interfacing with PEST(++)
BSD 3-Clause "New" or "Revised" License
169 stars 94 forks source link

Issues in parameterizing a MODFLOW6 timeseries file #424

Closed Paszka1 closed 10 months ago

Paszka1 commented 1 year ago

After trying a couple of options to parameterize a MODFLOW6 timeseries file I got stuck.

  1. The first attempt was to use PstFrom.add_parameters but it doesn't seem to work. The problem is with the columns as apparently the index_cols cannot be a list with a single integer. I have 1 (time) + 69 (well pumping rates) columns in the file and adding index_cols=[0] expects only one column to parameterize. Providing a list to use_cols doesn't fix the problem either.
  2. The second attempt was to create a template file using bits and pieces from the PstFrom library as I was hoping that the metadata will be taken into account by pst.add_parameters. The template file is built, parameters are added (though an error message is given; see below), but there are two issues then: (1) as noted, I added parameters with pst.add_parameters, and the multiplier info CSV file was not updated and obviously my hopes that the metadata will be considered flew away, and (2) I'm not sure what is the importance of the line in the template file starting with ,sidx,idx_strs,... just above the first row of the timeseries table in the template file. This is something pst.add_parameters would not understand giving the error message mentioned previously, so I wonder whether it could simply be omitted from the template file, or not?

Note, my model uses a DISV grid, but I don't think this has anything to do with the above. I don't know if it would make sense to add a add_parameters_from_tpl method to PstFrom?

Any thughts on this matter are greatly appreciated.

Paszka1 commented 1 year ago

Sorry guys.

I think I got the fix. Just needed to read again the documentation for PstFrom.add_parameters().

briochh commented 1 year ago

No problem @Paszka1. There are a lot of options in that add_parameters() call, it is easy to mis-step. Would it be possible to update this issue with your solution? This might be a fairly common use case, and I am sure you wont be the only one to hit a similar issue.

Paszka1 commented 1 year ago

Sorry guys.

I think I got the fix. Just needed to read again the documentation for PstFrom.add_parameters().

Paszka1 commented 1 year ago

Sure @briochh. The MODFLOW6 standard timeseries file has timeseries for 69 MAW wells, which means that I have 69 columns of monthly extraction rates plus a time column (column 0). What I first missed was that for parameter names, parameter groups and ultimate upper and lower bounds I needed to specify a list-like parameter with length equal to the length of use_cols. So the code section looks like below. The header of the TS file in my case has 7 rows (including a comment line added by Flopy) that need to be skipped. Parameters are by default defined as multipliers and log transformed.

I have a question here though. Although I'm not planning to use the ult_lbound and ult_ubound parameters, I'm not sure if my interpretation as seen in the code section below is correct. I.e., these parameters are only mentioned in the pyEMU documentation (not in the PEST++ documentation) and I understand that they are the actual parameters in the model input file following multiplication, which to me means that they can also be negative even for a log transformed parameter. Is this correct, or not? Please advise. The reason I'm asking is that the above interpretation seems to be conflicting with what is provided in the description of ult_lbound in the pyEMU doc.

Let me add that pyEMU is extremely useful and the GMDSI tutorials are also a great help. The documentation just needs to be read carefully.

use_rows = np.arange(sim.tdis.nper.array).tolist()

nwells = gwf.maw_1.nmawwells.array
use_cols = np.arange(1, nwells+1).tolist()

name = ["mawprategr"] * nwells
ult_lbound = [-500.0] * nwells # see my question above
ult_ubound = [-1.0] * nwells # see my question above

pf.add_parameters(ts_file,
                  par_type='grid',
                  geostruct=None, # it doesn't really make sense to add spatial or temporal correlation for extraction rates
                  par_name_base=name,
                  pargp=name,
                  mfile_skip=7,
                  index_cols=[0],
                  use_rows=use_rows,
                  use_cols=use_cols,
                  lower_bound=0.25,
                  upper_bound=2.5,
                  ult_lbound=ult_lbound,
                  ult_ubound=ult_ubound,
                 )