Simple PEST Calibration for MF6 model

[sebmeng]

Dear pyEMU community:
I wanted to try pyEMU in a simple steady state modflow 6 model to show the implementation to my advisor and continue using it in the transient calibration. The model has 4 layers, and in every layer there are zones defined for the hydraulic conductivity to be parameters for the PEST calibration. I also want to add the conductance of my river cells as calibration parameters. I have 32 points with observations in different layers to be added as targets for the PEST simulation. Here is my code:

import pyemu
import flopy
import flopy.utils
import numpy as np
import pandas as pd

#first i stored the observations as the format required by pyEMU:
# Select the columns and transpose the DataFrame
df_transposed = df_merged[['ID','Mean_phrea']].set_index('ID').transpose()

# Convert the column names to uppercase
df_transposed.columns = df_transposed.columns.str.upper()

# Add a new column "time" at the beginning of the DataFrame
df_transposed.insert(0, 'time', 1)

# Save the DataFrame to a CSV file without the index
df_transposed.to_csv('heads.csv', index=False)

df_transposed.head()


#Then I created zone_files for each layer:

# Get the number of layers
nlay = gwf.modelgrid.nlay

# Create an empty 3D array for the zones
zones = np.zeros((nlay, gwf.modelgrid.nrow, gwf.modelgrid.ncol), dtype=int)

# Initialize a dictionary to store the K values for each zone
zone_k_values = {}

# Loop through each layer and create zones
zone_counter = 1
for layer in range(k.shape[0]):
    unique_k_values = np.unique(k[layer])
    for k_value in unique_k_values:
        zones[layer][k[layer] == k_value] = zone_counter
        zone_k_values[zone_counter] = k_value  # Store the K value for the zone
        zone_counter += 1


#define the spatial reference and the PstFrom object:

sr = pyemu.helpers.SpatialReference.from_namfile(
        os.path.join( "Gondebeek_SS.nam"),
        delr=gwf.dis.delr.array, delc=gwf.dis.delc.array)

# Create a PstFrom object using pyEMU to generate the PEST control file and template files.
pf = pyemu.utils.PstFrom(original_d=".", new_d="pest_simulation", remove_existing=True,
                         longnames=True, spatial_reference=sr,
                         zero_based=False, start_datetime="1-1-2015")

#add the heads previously stored:

df = pd.read_csv("heads.csv",index_col=0)
df


# Add the observations to the PEST control file
hds_df = pf.add_observations("heads.csv", insfile="heads.csv.ins", index_cols="time",
                             use_cols=df.columns.values.tolist(), prefix="hds",)

hds_df

#then add the hk as parameters to the pf object:

# Add the zones to the PstFrom object
for i in range(1, zones.max() + 1):
    for layer in range(zones.shape[0]):
        zone_file = "zone_{0}_layer_{1}.dat".format(i, layer)
        zone_file_path = os.path.join("pest_simulation", zone_file)
        np.savetxt(zone_file_path, zones[layer] == i, fmt="%d")
        pf.add_parameters(filenames=zone_file, par_type="constant", par_name_base="zone_{0}_layer_{1}".format(i, layer),
                          pargp="zone_{0}_layer_{1}".format(i, layer))

And then i'm stuck :( I have some questions:

• How should I add the template file for the NPF package, where i want that PEST changes the values of the zones of HK?
• Were should I define the boundaries and initial values for my HK zones?
• How do I manipulate the observations? Using the groundwater utilitites (i have the files for using mf6mod2obs, but here the observations are handled in a different format than the neccesary for that utility)
• How do I add an instruction file for the observations?
• How do I run PEST through pyEMU? it is just for the creation of the pst simulation?

Maybe they are basic questions but any help would be really appreciated :)

If is needed I can share the files that I used for my simulation,

Thanks in advance!!!

Sebastian

[wkitlasten]

First, do you have all of your MODFLOW files written as external files? If not, I suggest that as your first step.

If so, adding a multiplier to k for each zone in each layer should be something like:

pf.add_parameters(filenames=exernal_k_filename, zone_array=zone_file, par_type='constant', par_name_base=your_str, pargp=your_str

You could also/alternatively add pilotpoints for each zone using par_type='pilotpoint' and passing pp_space and geostruct.

You can define the boundaries in the add_parameters() call with keywords: upper_bound=value, lower_bound=value

When you pfadd_observations() the file you point to needs to have the same format as you expect pest to read. So if you are not post-processing your mod2obs outputs at all, you need to point to the mod2obs output file with pf.add_observations(). Then you will need to assign the actual (real world) value you want pest to history match to. The .ins file is created with the pf.add_observations() call.

Or, you can build your instruction files yourself and use: add_observations_from_ins()

You can run pest from pyemu with pyemu.helpers.start_workers()

Hopefully that gets you started. Check the docs/help for heaps of additional tips: help(pf.add_parameters), etc