PEST parameters work by modifying the values assigned to data sets. For example if a data set is assigned a value of 2 when parameters are not used, then if a PEST parameter with a value of 3 is used with the data set, the final value for the data set will be 2*3 = 6.
When PEST is calibrating a model, it does so by changing some of the model inputs (parameters). As described in detail in Doherty and Hunt (2010), the decisions on which inputs become parameters influences model calibration and how the model can be used (e.g., uncertainty analysis). The user decides which model inputs PEST can vary by first defining parameters in the Manage Parameters dialog box or for some MODFLOW packages in the MODFLOW Packages and Programs dialog box and then using the parameters in the Data Sets dialog box or the Object Properties dialog box.
For MODFLOW-2005 and related models, parameters can be defined for some packages using capabilities built directly into MODFLOW. Built-in parameters have been eliminated from MODFLOW 6 but ModelMuse still allows you to use MODFLOW-2005 style parameters for some packages by generating MODFLOW 6 input files based on how parameters were used with MODFLOW-2005. For SUTRA, MODFLOW 6, MODFLOW-2005 and related models, the user can also define parameters that are unrelated to any package. Their parameter type is identified as "PEST" in the Manage Parameters dialog box.
Both MODFLOW-2005 style parameters and PEST parameters can be used with PEST. However, the way that PEST parameters are applied to data sets is somewhat different from how MODFLOW-2005 style parameters are applied. With MODFLOW-2005 style parameters for array data, the user can define zone and multiplier arrays. If a zone array is used, the parameter will be applied to cells where the zone array data set is set to true. If a multiplier array is used, the value applied to a cell is the parameter value times the multiplier array value. It is possible to have more than one MODFLOW-2005 parameter apply to the same cell of the same data set by setting the zone arrays for two or more parameters to true for the same cell. In that case, the value applied to the cell would be the sum of the values that would be applied by each parameter individually.
If PEST parameters are to be assigned to a data set, the “PEST Parameters used” checkbox is checked on the PEST Parameters tab of the Data Sets dialog box. The PEST parameters tab will only be present for some data sets. If it is not present, PEST parameters cannot be used with that data set. ModelMuse only allows PEST parameters to be used with real-number data sets. It does not allow PEST parameters to be used with data sets that define the layer structure of the model despite those being real-number data sets. In addition, in SUTRA models, parameters can only be applied to data sets corresponding to PMID, PMIN, ALMID, ALMIN, ATMID, and ATMIN if anisotropy has not been selected for those data sets on the Anisotropy pane of the SUTRA Options dialog box.
For every data set for which the “PEST Parameters used” checkbox is checked, a corresponding text data set will be created when the Apply button in the Data Sets dialog box is clicked. The name of the corresponding data set will be the same as its parent data set with "_Parameter_Names" appended. This "_Parameter_Names" data set is used to designate the locations where a PEST parameter is applied. A PEST parameter will be applied to any cell in the parent data set for which the value of the _Parameter_Names data is set to the name of one of the PEST parameters. For example, suppose in a MODFLOW model, there is a PEST parameter named “MyKx” and the “PEST Parameters used” checkbox is checked for the Kx data set. Then a polygon object is used to set the value of the Kx_Parameter_Names data set to “MyKx” for some cells. When the model was run the MyKx parameter would be used with the Kx data set in the cells for which Kx_Parameter_Names equaled “MyKx”. At other cells, different PEST parameters or no PEST parameter at all might be applied. The Formula Editor dialog box lists all the PEST parameters, and these can be used as the formula for a _Parameter_Names data set either as the default formula in the Data Sets dialog box or in the Data Sets tab of the Object Properties dialog box. When a PEST parameter is applied to a cell or element in a data set, the value at that cell or element is multiplied by the parameter value.
Using Pilot Points
A PEST parameter can be associated with pilot points (Doherty, 2003; Doherty and others, 2010). Pilot point approaches used in ModelMuse are explained in the manuals for PLPROC and the PEST Groundwater Utilities available from the PEST home page (https://pesthomepage.org/). ModelMuse uses PLPROC to implement pilot points. Section 4 of the documentation for the PEST Groundwater Utilities provides a conceptual description of pilot points and how they are used. In brief, pilot points are a parameterization device that facilitates increased parameter flexibility by estimating properties at user-specified locations in the grid; the remainder of the grid is then filled by kriging. Pilot points can be grouped with zonation, but such zonation is not required. For instance, suppose the hydraulic conductivity of an aquifer is being estimated. From aquifer tests, it may be known that the hydraulic conductivity varies from place to place but the spatial variation might not be well defined. One way to approach this would be to specify zones within the aquifer and to have PEST estimate separate uniform parameter values within each zone. A potential drawback of this approach is that the chosen zonation might not be optimal. Pilot points provide a way to get around this problem. The modeler designates a number of points and assigns a value to each point. Between points, values are assigned to each cell via kriging. Within PEST, the value assigned to each pilot point is a parameter that can be adjusted by PEST to improve the fit between the observed and simulated values. If all pilot points are tied in the PEST control file, they will act like a zone. Likewise, a zone with zero or on associated pilot points will also produce zone-like results. Such aspects can be useful in stepwise modeling, where model complexity is added in sequential steps in response to model performance.
Each pilot point value is treated as a separate parameter in PEST. PEST will typically run the model at least once for each parameter during each parameter estimation. Thus, the number of pilot points can have a large effect on the time required for parameter estimation. You can estimate the time required for one parameter estimation by multiplying the number of parameters by the time required to run the model one time. You can find the time required to run the model by determining the time difference between when the last model input file was created by ModelMuse and the time the last output file was created by the model. The last input file for the model created by ModelMuse is typically "RunModel.bat". The number of parameters is the variable NPAR in the PEST control file. NPAR is the first number outside of a comment in the control data section of the PEST control file.
Page 9 of Doherty and Hunt (2010) provide suggestions for pilot point placement. Candidate pilot point locations are defined in the Pilot Points pane of the PEST Properties dialog box. They can be specified in several ways. They can occur in a regular pattern, be between point observations, be specified individually, or be any combination of these methods. The modeler also specifies a pilot point buffer that affects how pilot points are defined.
Whether or not a particular parameter is associated with pilot points is specified in the Manage Parameters dialog box. If pilot points are used with a parameter, the treatment of the parameter is somewhat modified while the parameter estimation process is running. Instead of multiplying the data set value by the parameter value, each pilot point will be an independent parameter. The PEST utility program PLPROC is used to interpolate among the pilot point values, and the interpolated values will be assigned to the data set wherever the parent parameter is used. The pilot points used with a particular parameter for a particular data set will be any of the candidate pilot points defined on the Pilot Points pane of the PEST Properties dialog box that are no farther away than the pilot point buffer from a cell center (MODFLOW) or node location or element center (SUTRA) that is part of the zone where the parameter applies. The initial value assigned to a pilot point will be the value in the corresponding data set in ModelMuse. However, if the parent parameter for a pilot point is not used for the pilot point location, the nearest location for which it is used will supply the initial value for the pilot point.
To clarify how ModelMuse attributes initial values, consider a simple model with 10 rows, 10 columns, and 1 layer. The rows and columns have a spacing of 100 meters. The model has two parameters defined - KLeft and KRight as illustrated below. They will be used to define the hydraulic conductivity in the X direction on the left and right halves of the model respectively. Both parameters have initial values of 1. Pilot points are used with KRight but not with KLeft.
Screen capture of the Manage Parameters dialog box with two parameters defined.
There are 25 candidate pilot points defined that are spaced 200 meters apart as illustrated below. Nine of the pilot points are outside the model grid.
Screen capture illustrating the locations of the candidate pilot points.
The pilot point buffer in this model is 290 m. Because the leftmost column of pilot points is more than 290 meters from the center of any cell in the right half of the model (where pilot points are applied) that column of pilot points will not be used. On the other hand, the pilot points above and to the right of the grid (except the one in the first column) are within the pilot point buffer and will be used. Their initial values will come from the cell to which each of them is closest. In addition, the second column of candidate pilot points will be used because they are within 290 m of the right-hand side of the model where the KRight parameter is applied. Their initial values will come from the nearest cell on the left-hand side of the model.
The model has specified heads in the first and last columns (0.1m in column 1 and 1.0 m in column 10). The model has four head observations: 0.45 m in row 2, column 3; 0.75 m in row 2, column 6; 0.35 m in row 8, column 3; and 0.65 m in row 8, column 6.
After PEST has finished running, the values assigned to the Kx data can be displayed using "File|Import|Gridded Data Files..." To do so, select the file containing the final values in the "arrays" subdirectory of the directory in which the model ran. The file will have the extension ".arrays." The final distribution of Kx is shown below. Note that the left half of the model where no pilot points were used has a uniform distribution of Kx values whereas in the right half where pilot points were used, the hydraulic conductivity varies among the cells.
Diagram showing distribution of Kx values after PEST estimated parameters.
It is recommended that when pilot points are used either reasonable values should be assigned to inactive cells for all estimated data sets or else the parameter name should be assigned a value of "" in all inactive cells. For example, suppose Kx is being estimated, Kx is assigned a value of zero in inactive cells, and pilot points are defined outside the grid or in an inactive cell. For the pilot points in inactive cells and perhaps some of the cells outside the grid, the initial value will be zero. The lower and upper bounds for those pilot points will also be set to zero. The upper and lower bounds for the pilot points will be set as the initial value times the ratio of the upper or lower bound of the base parameter times the initial value of the base parameter. That will result in the upper and lower bounds for those pilot points also being set to zero.
Important Caveat
The formulas assigning values to data sets are only applied in ModelMuse. When PEST is running, those formulas are not used. This is especially important to keep in mind when a data set is related to a direction such as the Kx, Ky, and Kz data sets in MODFLOW.
PEST is a universal parameter estimation code that can be applied to most any forward run model. It achieves this wide application because it operates – that is, adjusts parameters and evaluates outputs – outside of the forward model code itself. Therefore, formulas assigning values to data sets are only applied in ModelMuse to create the files run by PEST. This is especially important to keep in mind when a data set is related to a direction such as the Kx, Ky, and Kz data sets in MODFLOW. Having PEST assign values to the Kx data set does not mean that the Ky and Kz data sets will be automatically updated to the MODFLOW forward run too. Moreover, when Kx, Ky, and Kz are specified as independent parameters there can be no expectation that geologically reasonable anisotropy is maintained as the parameter estimation proceeds. As far as PEST is concerned, those are independent properties unless specifically linked through parameter pre-processing outside of ModelMuse (e.g., PEST utility PAR2PAR – see PEST manual for additional discussion). For MODFLOW models, there are options to estimate horizontal and vertical anisotropy instead of independently estimating Ky and Kz – this allows the use of parameter bounds to enforce geologically realistic anisotropy during the exploration of parameter space. For MODFLOW 6, these are options in the Node Property Flow package (K22OVERK and K33OVERK). For the LPF and UPW packages, horizontal anisotropy is used automatically but vertical anisotropy can be specified by an option on the Basics tab of the MODFLOW Layer Groups dialog box. For the BCF package, horizontal anisotropy is used automatically. Vertical leakance in the BCF package is a function of the vertical hydraulic conductivities of the more than one layer. In the Hydrogeologic Unit Flow (HUF) package, all data are specified using parameters. You can define parameters for horizontal and vertical anisotropy in the HUF package. When working with the HUF package, the values assigned to cells can be a composite value from several different hydrogeologic units. Values for individual hydrogeological units can be generated using HUFPrint (Banta and Provost, 2008) or functions built into ModelMuse for comparison with the conceptual model.
With SUTRA models, having PEST assign values to the dispersivity, permeability, or hydraulic conductivity in the "max" direction does not mean PEST will assign values in the "mid" or "min" directions unless anisotropy is used. ModelMuse allows you to use horizontal and vertical anisotropy for those data sets and anisotropy is used by default. The anisotropy options are specified on the Anisotropy pane of the SUTRA Options dialog box. When these options are used, ModelMuse generates a template file for the main SUTRA input file that uses a formula to relate the "mid" or "min" data sets to the corresponding “max” data set.
