Navigation:  Quick Start Tutorials >

PEST with Pilot Points Tutorial

 

This section provides instructions on using Visual MODFLOW Flex to setup, run, and interpret a Parameter Estimation/Predictive Analysis simulation. In addition, this tutorial provides a brief description of the input parameters and settings required by PEST. A detailed description of the algorithms, parameters, input files, and other options for PEST are available in the PEST User Documentation, which are available at the PEST website (www.PESTHomepage.org). 

This exercise demonstrates some of the advanced and exciting opportunities for model calibration and uncertainty analysis using PEST in Visual MODFLOW Flex. This exercise is based on the problem described in "Using Pilot Points to Calibrate a MODFLOW/MT3D Model" by John Doherty (2008), and has been adjusted to work with the PEST workflow inside Visual MODFLOW Flex.

 

Before you start

You are encouraged to familiarize yourself with the concepts and applications of PEST prior to using in Visual MODFLOW Flex. The time spent on this will make your experience with parameter estimation much more productive, and will likely help you to overcome any difficulties you may experience the first time you run PEST.

 

In addition, if you are not familiar with the Visual MODFLOW Flex graphical environment, please take a few minutes to review the Program Overview section.

 

Note:

You must have a Pro or Premium license in order to use the PEST module in Visual MODFLOW Flex.

 

Objectives

By the time you have finished this exercise you will have:

used pilot points as a means of characterizing the spatial distribution of an aquifer hydraulic property;

used PEST’s advanced regularization functionality in conjunction with geostatistically-based regularization constraints;

used the Visual MODFLOW Flex GUI to:

oBuild input files for PEST

oRun PEST

oAnalyze the results

oSave adjusted parameters as new model inputs

 

Required Files

Several files are required for this exercise, which are included with the Visual MODFLOW Flex installation. These files are available in the installation documents folder, typically:

       C:\Users\Public\Public Documents\Visual MODFLOW Flex\Tutorials\PEST

 

If you cannot find these files, please download the PEST Tutorial project from our website, and unzip to a desired folder on your computer.

 

Introduction

The groundwater model that you will use in this tutorial has already been built. In the first part of the tutorial, you will open the project, translate and run the model, and start a new PEST workflow:

Launch Visual MODFLOW Flex

Click File > Open Project

Navigate to your Public 'My Documents' folder, and locate '\Documents\Visual MODFLOW Flex\Tutorials\PEST', and open the 'pest-tutorial.amd' file

Open, to open the project

In the Numerical Workflow window, locate the workflow tree (this is shown on the left in the window below)

If you're not already on this step, click on the Select Run Type step in the numerical workflow tree.

Select the Single Run button

 

 

Under the Flow Engine dropdown menu, select MODFLOW-2005 and select the Next Step button []. You will now be on the Translate step.

Translate and Run the model. Once the model converges, return to the Select Run Type step. Note that a model must be ran at least once before PEST can be run.

 

Please Note: When using PEST, the head close change criterion (typically called HCLOSE or similar in the solver settings) should generally be as small as practical to increase the precision of the derivative calculations. For this tutorial we will leave them on their defaults to decrease model run times.

 

Select the PEST Run button to open the PEST workflow window.

A new PEST Workflow window will load as shown below:

 

Please Note: if this does not appear right away, select this tab from the list of workflows across the top of the window, highlighted in image below)

 

 

Define Observations

The first step in the PEST workflow is to define the observations you want to include for calculating the objective function, and assign weights to various observations. In this example, there are 21 observation wells where heads were measured at several intervals.  These observations were assigned in the main modeling workflow (i.e. the "Sample2grid-Run" workflow) and are automatically added in the PEST workflow. All head observations are selected by default, and will be used with a default weight = 1.

No changes will be needed at this step for this exercise

Click [] (Next Step) to proceed to the 'Define Property Parameters' step.

 

Please Note: if you receive a pop-up warning that "Observation settings have been changed...", simply select OK

 

Define Property Parameters

At the Define Property Parameters step in the workflow, you will define which parameters will be included in the objective function, that is which parameters will be adjusted by the PEST algorithm to reduce differences between observed head values and corresponding simulated head values.

The 'Define Property Parameters' window will appear:

 

 

The table at the top of this step allows you to select which input model properties and boundary conditions you want to include as parameters in the objective function. For each parameter, you can specify whether to "Tie" it to another parameter so that it will be varied during the PEST run in relation to the selected anchor parameter (for example, if Kz is tied to Kx, it will be varied in the PEST run so that the ratio between Kx and Kz at each pilot point remains constant). You can also specify a Transformation option (by default, all Conductivity parameters are set to Log transformation) for the parameter sensitivity derivatives.

The table at the bottom of this step allows you to select which model parameters you want to vary as part of the PEST run. For this exercise, all Kx property zones will be included. You can also specify bounds on the values for each parameter zone:

Check Kx on the upper table, to include it as a PEST parameter

Enter 1 for the Minimum for each zone in the lower table

Enter 300 for the Maximum for each zone in the lower table

 

 

Click (Next Step) to proceed to the Define Pilot Points step.

 

Define Pilot Points

At this step, you will define pilot points, which are the point locations in the model where PEST will estimate the parameters and spatially interpolate them. One set of pilot points is required for each distinct parameter zone: a single pilot point in a parameter zone will result in a homogeneous parameter field for that zone and multiple pilot points in a parameter zone will result in a heterogeneous, continuously distributed parameter field for that zone.  The pilot point approach implemented in PEST allows you to methodically determine sets of parameters for your groundwater model that minimize the PEST objective function (e.g. reduce the global differences between observed and simulated values).

The Define Pilot Points window will appear:

 

 

Pilot points are XY points with an initial value for each parameter. Pilot points can be imported from .TXT file, XLS, .SHP file, or assigned manually by digitizing in the 2D environment. An example of the pilot points on top of conductivity zones is shown below:

 

 

The general steps to assign pilot points in the PEST workflow are as follows:

 

Select pilot-points-zone1 data object from the Data Explorer (tree)

Click [] (Insert button) at the top of the 'Define Pilot Points' window to add these points.

Select which parameter zones the points represent under the Parameter Zones grid in the upper right section of the window.

For this set of pilot points, Kx-Zone1 is fine (this should be selected by default)

Your display should now appear as shown below:

 

 

Repeat these steps to define the pilot points for Kx Zone2, Kx Zone3, and Kx Zone 4:

Select 'pilot-points-zone2' data object from the Data Explorer (tree)

Click [] (Insert button) at the top of the 'Define Pilot Points' window to add these points.

Select Kx-Zone2 for these pilot points (if it is not already selected)

 

Select 'pilot-points-zone3' data object from the Data Explorer (tree)

Click [] (Insert button) at the top of the 'Define Pilot Points' window to add these points.

Select KxZone3 for these pilot points

 

Select 'pilot-points-zone4' data object from the Data Explorer (tree)

Click [] (Insert button) at the top of the 'Define Pilot Points' window to add these points.

Select Kx-Zone4 for these pilot points

 

When you are finished, your display should appear similar to the figure below.

 

 

In the table at the bottom of the window, you can adjust some parameters for specific pilot points. Above this table, there is a combo box that allows you to select which set of pilot points is shown in the table based on the selected parameter zone.

The main use of this table is to specify initial values for the pilot points in each parameter zone and which of these pilot points will be used and/or set as "Fixed" during subsequent PEST runs. Fixed pilot points are those locations where you are confident in the measured parameter value (e.g. based on pumping/slug test locations) and you want these values to remain unchanged during the PEST run.

In this example, leave the defaults as is.

Click [] (Next Step) to proceed to the Define Kriging Variograms.

 

Define Kriging Parameters

The use of pilot points in characterizing the spatial distribution of a hydraulic property must be accompanied by a mechanism whereby hydraulic property values assigned to pilot points are spatially interpolated to the cells of the finite difference grid. In PEST, this spatial interpolation is accomplished using the Kriging algorithm. Kriging is a method of spatial interpolation from the field of geostatistics. The cornerstone of geostatistics is the variogram. A variogram describes the extent to which hydraulic property values (or any other type of data) pertaining to any two points are likely to be different from each other as a function of the distance between those points.

One of the benefits of using Kriging as a basis for spatial interpolation is that the factors by which hydraulic properties at pilot points are multiplied before summation to obtain the hydraulic property value at a particular grid cell are independent of the actual hydraulic property values at the pilot points. Hence a set of “Kriging factors” pertaining to each of the cells of the finite difference grid can be calculated in advance of the actual interpolation process. Since the interpolation process is undertaken again and again as the model is run repeatedly by PEST, the fact that it is not necessary to repeat the calculation of the Kriging factors on each occasion that the model is run can result in large savings in the time required to complete the overall parameter estimation process.

In this example, there are 4 hydraulic property zones as shown below:

 

 

Zone 1: River Alluvium

Zone 2: Creek Alluvium

Zone 3: Western Basalt

Zone 4: Eastern Basalt

 

Each zone is represented by a geostatistical structure. Each of these structures cites one variogram (though it could cite up to five, one for each groundwater flow property: Kx, Ky, Kz, Ss, Sy). In this tutorial:

“Structure1” will be used to characterize Zone 1 of our model domain (i.e. the river alluvium),

“Structure2” will be used to characterize Zone 2 (i.e. the creek alluvium),

“Structure3” will be used to characterize Zone 3,(i.e. the western basalt), and

"Structure4” will be used to characterize Zone 4 (i.e. the eastern basalt).
 
Note that the variogram assigned to these zones 3 and 4 is not important: because there is only one pilot point assigned to each of them, all cells within these zones will be assigned the one interpolated value (same as the respective pilot point) irrespective of the variogram. Assigning a single pilot point to a given zone is effectively the same as the approach that was standard in older versions of PEST where hydrogeologic properties were homogenous within a zone.

 

In Visual MODFLOW Flex, a default variogram is generated for each parameter zone, based on a type 2 (Exponential) distribution.

 

 

However, the structures should be verified/modified.

Uncheck the Calculate variograms automatically checkbox

Click 'Zone1' under 'Kx' in the tree. The display will appear as below:

 

 

In the editor on the right side of the display, for 'Transform', set this to log, if it isn't selected already.

Repeat these steps, setting log 'Transform' for Kx in Zones 2, 3, and 4.

 

Any variogram cited in each of these structures will pertain to the spatial distribution of the logarithm of the pertinent hydraulic property. This is in accordance with the fact that most studies cited in the groundwater literature which treat transmissivity and/or hydraulic conductivity as a regionalized variable indicate that its distribution is better described by a log variogram than by a variogram based on native property values.

For this example, the default parameters for the variograms (all will use Exponential) is sufficient. However, for Zone2, we will define a value of 2.0 for the Anisotropy, with the direction of anisotropy coinciding with the direction of the creek. Alignment in the direction of the creek is based on the premise that channel structures within this old creek valley will make it more likely for hydraulic property similarity to prevail in this direction than in a direction at right angles to it.

The variogram parameters can be adjusted by selecting a Variogram from the tree as shown below:

Click 'Variogramkx2' under the Variograms node in the tree.

Enter 2 for the 'Anisotropy'.

Enter 45 for the 'Bearing'. This value allows you to make the variogram anisotropic in a certain direction; 'Bearing' is an angle of rotation specified in degrees counter-clockwise from the positive X-Axis (East).

The display should appear as below:

 

 

Click [] (Next Step) to proceed to the Select PEST Variant step

 

Select PEST Variant

In this step, you have the option of choosing which PEST variant to use. Currently, the PEST engines PEST, PEST-HP, and PEST-HP (MKL) are supported by Visual MODFLOW Flex. See the PEST variant documentation page for a detailed discussion of the differences between each. Importantly, the basic PEST engine runs models as a serial process (i.e. sequentially, one at a time), whereas PEST-HP and PEST-HP (MKL) can run multiple models in parallel. PEST-HP can dramatically decrease the length of a PEST run, especially for computers with many cores.

That said, for this tutorial, we will use the basic PEST engine, which should be selected by default. You may want to experiment with using the PEST-HP engines to see how they compare.

Select PEST from the PEST Variant dropdown menu, if not selected already

Click [] (Next Step) to proceed to the PEST Control Settings step

 

PEST Control Settings

This step contains settings for the PEST control file. In the PEST Controls Settings section, the options for the Control Data section in the PEST control file are available. These options are separated according to which line they correspond to in the block; a full description of each variable may be found in Chapter 4 of the PEST manual (Doherty 2023a). Among other things, these variables affect how much PEST can change variables between iterations, the target value of the objective function phi, and how PEST handles model run failures. These settings may have a large impact on the speed and stability of the PEST run, as well as the possibility that PEST produces over-fitted parameter values.

The Singular Value Decomposition (SVD) section contains options for singular value decomposition regularization, which can improve the stability of working with ill-posed problems with non-unique solutions, as is often the case with groundwater models which typically have more uncertainty associated with inputs (i.e. parameters) than with outputs (i.e. observations). Note that this is distinct from "SVD Assist" regularization, which is discussed below. SVD achieves numerical stability in PEST by subtracting parameters (or combinations of parameters) from the calibration process. As a result of the subtraction, the calibration process is no longer required to estimate values for those individual parameters (or combinations of correlated parameters) that are inestimable based on the contribution of their values to improving the objective function (e.g. fitting the calibration dataset). These combinations are automatically determined through singular value decomposition (SVD) of the weighted Jacobian matrix. As before, a detailed description can be found in chapter 4 of the PEST manual.

 

We will leave the settings on their defaults as they work well for this model, and proceed to the next step.

Click [] (Next Step) to proceed to the Select Run Type step

Select Run Type

At this step, choose the type of PEST Run. If you want to run PEST, then some additional options will be needed, such as define regularization and adjust the PEST control file.  If you want to run Sensitivity Analysis, this can also be launched.

 

Click the [Parameter Estimation] button from the main window, as shown above.

The next step will be to choose the Regularization options.

 

Select Regularization

At this step, you can choose the type of Regularization to run.  Regularization is a process whereby additional information is introduced into the objective function of the model in the form of ‘prior information equations’. These prior information equations can take many forms:

Tikhonov Regularization can induce PEST to prefer parameter values which are as close to their initial values as possible. Thus, Tikhonov regularization method adds more information to the overall problem, but this additional information helps to guide PEST toward acceptable parameter estimates as informed by your knowledge of the modeled area, helps to avoid ‘over-fitting’ solutions, and reduces the non-uniqueness of potential solutions.
 

Singular Value Decomposition (SVD) Assist achieves numerical stability in PEST by systematically reducing the parameter space initially composed of independent parameters with their own dimension to a lower dimensional space of orthogonal super-parameters (where each super-parameter is composed of parameters or combinations of parameters).  This reduction in parameter space significantly lowers the computational burden of running PEST, since the core of parameter estimation evaluations is calculating a Jacobian matrix which requires at least one model run per parameter (in the case of a standard PEST run) or at least once per super-parameter (in the case of SVD-Assist). The combinations of parameters that make up each super-parameter are automatically determined through singular value decomposition (SVD) of the weighted Jacobian matrix, which is calculated by performing a Sensitivity Analysis as a prerequisite.

 

For this exercise we will not perform any regularization.

 

Click No Regularization button from the main window

The next step will be to adjust the PEST control file.

 

Edit PEST Run Settings

The last step before running PEST is to review and adjust the PEST Control file.

 

 

If you are familiar with the PEST file format/structure, you can adjust the PEST Control file in this window, or copy into a text editor, make changes, and paste the adjusted contents back in this window.  A full explanation of the PEST control file is available in the PEST manual. For this exercise, the default values are fine.

Before starting the PEST Run, it is a good idea to check the PEST Input files.  PEST provides a utility to do this, called PESTCheck.

 

Click the [] button on the workflow toolbar

You should receive a confirmation that no errors were found.

Click [OK] to proceed.

Click [] (Next Step) to proceed to the Parameter Estimation.

 

Parameter Estimation

At this step in the workflow, you will run PEST:

Click [] (Run PEST) to start the PEST Run.

 

The PEST program (pest.exe) will load in a DOS command window, and show the progress as seen below.

 

The PEST run should several minutes. For unregularized models, the number of model runs required for each iteration will be strictly greater than the number of parameters, which is why parallel model runs with PEST-HP can greatly reduce the amount of time required for a full PEST run. The groundwater model used in this tutorial is relatively small and has relatively few parameters - PEST runs involving larger models with longer run times and many parameters may take considerably more time: up to several days or even longer.

As PEST runs, you should see the objective function (phi) decreasing over each optimization iteration; pay attention to these values in the DOS window. In this tutorial, PEST will run a total of about 20 optimization iterations and a total of approximately 1200 model runs. PEST should reach a final objective function (phi) value of approximately 1E-02.

When PEST finishes, you should see a confirmation message in the main window, below the 'PEST Run Log' tab, as shown below.

 

After the PEST run completes, you can analyze the results.

Click [] (Next Step) to proceed to the Analyze Results step.

 

Analyze Results

Visual MODFLOW Flex will present the results of the PEST run, with one tab per output file:

 

 

Record file (.REC): contains parameter values, objective function, sensitivities, etc..

Sensitivities for Observations (.SEO): contains the observed and simulated values with the sensitivities

Sensitivities for Parameters (.SEN): provides the information on parameter sensitivities

Residuals (.RES): contains the adjusted calculated vs. observed values and residuals

 

The results from these files can be exported into Excel for charting.

Click [] to export the PEST results to an Excel spreadsheet.

 

A window will appear which allows you to select the location of the exported Excel spreadsheet, and to select from available templates. Two templates are available; one which will export the available PEST data to a series of worksheets/tables, and a 2nd template which will also generate a number of charts for each type of PEST data.

Select the 'pest-results-chart-and-table.xlsx' template

Click the checkbox besides 'Open file when finished'

Select an output location

Click Export; the Export Pest Results window should look like the image below before exporting:

 

 

An Excel spreadsheet will open which summarizes the results of the PEST run. Several data tables and charts are included; spend a few minutes reviewing them. The screenshots below display a small selection of the available PEST data (yours may look slightly different for this model). For more information about the PEST results please visit the PEST website (www.PESTHomepage.org).

 

If the results look reasonable, you can save the adjusted Kx parameter zonation as inputs for a new model, as explained in the next section.

Click [] (Next Step) to proceed.

 

Save PEST Parameters as New Inputs

After reviewing the PEST output, if the adjusted parameter values seem reasonable, you can save these parameters as inputs for a new model run.

 

Click on the [Create New Model Run with PEST Results] button.

 

Visual MODFLOW Flex will save the adjusted model parameters in a new model run within the same project. This new model run will appear in the Model Explorer below the most recent model run. A new workflow window will also appear with this model run.

 

Click on Define Properties from the workflow tree

From the Toolbox, select Kx as shown below; you should then see a color flood of the Kx values.

.

 

You can mouse over the 2D display to see the range of Kx values. Click on the [] button on the toolbar to show color shading with contour lines.

You must Translate and Run this new model run in order to see the updated MODFLOW results using the adjusted Kx parameters from PEST.

 

*****This concludes the Model Calibration Using PEST with Pilot Points tutorial*****

 

 


Page url:http://www.waterloohydrogeologic.com/help/vmod-flex/index.html?vm_pest_tutorial.htm