Using PEST with AnAqSim

Overview of the Process

AnAqSim, an easy to use Windows-based analytic element groundwater flow model, can also be run from the command line or a batch file. This ability allows the user to move beyond AnAqSim’s built-in manual model calibration features (maintaining a calibration target data set, and plotting a
map of residuals) to employ the power of PEST for automated model calibration.

PEST is a general-purpose parameter estimation software tool that has been widely used for calibrating MODFLOW models and other groundwater flow models. To use PEST with AnAqSim the user is required to create a batch file to run AnAqSim from the command line and create a template file for the AnAqSim model and an instruction file to read the AnAqSim output into PEST.

The creation of a simple AnAqSim script file described in section 4 of the AnAqSim User Guide (http://www.fittsgeosolutions.com/AnAqSimUserGuide.pdf) allows the user to automate AnAqSim runs from the command line. Specifically the user can specify where the run log should be written and what analysis commands should be executed. A list of the analysis commands is included in
Section 2 below. Using this command line functionality for AnAqSim, an automated PEST calibration can be performed for any AnAqSim model. A typical PEST run will include the following user created files:

• The PEST control file (*.pst)
• Any template files (*.tpl) that PEST needs (for AnAqSim there is only one template file because all of the AnAqSim inputs are stored in the *.anaq file in xml format).
• Any instruction files (*.ins) for reading the output (for AnAqSim there is only one instruction file because all of the output will be written to the same *.out file).
•  (Optional) A parallel run management file (.rmf) if you wish to use multiple processors on your own computer or across a network of computers.

To get started, PEST  is freely available from the USGS website at https://www.usgs.gov/software/pest-parameter-estimation-code-optimized-large-environmental-models-version-3

AnAqSim analysis commands

• initialheadsfile <name.hds> (Optional) name of the initial heads file if the model is
  set up to use an initial heads file
• outputfile <name.out> file where AnAqSim will write the results of all
  subsequent commands
• headspecifiedwells writes the discharges of any head-specified wells to the output file
• dischargespecifiedwells writes the heads at discharge-specified wells to the output file
• headspecifiedlines writes the discharges of internal head-specified line boundaries to the output
  file
• rivers writes the discharges of river line boundaries to the output file
• calibration writes the calibration results to the output file (includes the modeled
  head at each calibration target)
• verticalleakage writes the vertical leakages over polygon areas to the output file
• headoutputfile <name.csv> writes the heads to a csv file based on the current contour settings and plot window
• exit instructs AnAqSim to close the current model window

Simple Example

To illustrate the use of PEST with AnAqSim we will present a simple model parameter fitting example. The following description is intended only to show the kinds of things that AnAqSim is capable of when coupled with PEST. It is not intended to be a step-by-step tutorial or set of instructions.

The example provided below is a simple box model that only requires the free educational version of AnAqSim, available from http://www.fittsgeosolutions.com/downloads-licensing/ (be sure to download the AnAqSimEDU version unless you have purchased a license for the full version of AnAqSim).

The box model is a 5000 ft x 5000 ft square area of homogeneous sandy aquifer material with flow from left to right and three pumping wells located near the center. The box model was initially set up with a specified hydraulic conductivity and pumping rates for the wells to calculate the hydraulic head field within the model. From within that head field, we recorded the head elevations at six observation well locations.

The goal of this example will be to have PEST work from the six recorded head values as its inputs to find the “unknown” pumping rates of the three wells and hydraulic conductivity of the aquifer. The model domain and observation well locations are shown in the figure below.

Starting Model

Starting AnAqSim Model

Observation wells (calibration targets) are added to an AnAqSim model by selecting Analysis Input/Calibration Targets/Head and entering the names and locations of the monitoring wells as well as the observed (measured) head at each location. For PEST runs, the observed head value entered is not important, only the locations of the observation wells. This is so we can extract the modeled heads from AnAqSim after each model run using the command line; the observed heads are input directly into the PEST .pst file and are not extracted from the AnAqSim model.

The first PEST file that we need to create is the PEST control file (ExtracWells.pst). A brief discussion of the relevant sections of the file will be included here, but we suggest that those interested in a more detailed description of the PEST input variables read the information contained
in the USGS report describing PEST++ Version 3 (https://pubs.usgs.gov/tm/07/c12//tm7c12.pdf), specifically Appendix 1.

PEST Files

The ExtracWells.pst file is shown below, select parameters that need to be changed to adapt this example to other AnAqSim models are emphasized in blue. The remaining parameters affect how the
calibration is carried out by PEST and are considered beyond the scope of this example.

ExtracWells.pst

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

pcf * control data restart estimation 4 6 1 0 1 1 1 single point 1 0 0 10.0 -3.0 0.3 0.03 10 10.0 10.0 0.001 0.1 100 0.005 4 4 0.005 4 1 1 1 * parameter groups kh relative 0.01 0.0 switch 2.0 parabolic * parameter data kh_1 log factor 100.0 1.0 300.0 kh 1.0 0.0 1 p_well1 none factor -4000.0 -25000.0 -1.0 kh 1.0 0.0 1 p_well2 none factor -4000.0 -25000.0 -1.0 kh 1.0 0.0 1 p_well3 none factor -4000.0 -25000.0 -1.0 kh 1.0 0.0 1 * observation groups heads * observation data cal1 89.34256 1.0 heads cal2 80.57272 1.0 heads cal3 97.98874 1.0 heads cal4 82.09399 1.0 heads cal5 93.70227 1.0 heads cal6 95.60841 1.0 heads * model command line run_anaq.bat * model input/output ExtracWells.tpl ExtracWells.anaq ExtracWells.ins run1.out

The first line, pcf, indicates that this is the “PEST control file”. The next section contains the control data variables. The first blue variable, 4, is the number of parameters that PEST will be trying to calibrate, the second blue variable, 6, is the number of observation points that PEST will use to assess each run’s performance. These two variables must agree with the number of entries for 

* parameter data and * observation data.

The entries in each row under * parameter data contain the unique name of each parameter that will be adjusted during calibration; whether the variable should be adjusted using a log transform or not (for hydraulic conductivities, where values can range over several orders of magnitude, a log transform is recommended); the starting value that the parameter will take; the allowable minimum parameter value; and the allowable maximum parameter value.

The entries under * observation data contain the unique name for each observation point; the measured value for each observation point; the weight of each observation point; and what kind of observation it is (this can be named anything).

The command beneath * model command line is the command or batch file that is run for each model call. In this case, we are using the batch file run_anaq.bat.

The commands beneath * model input/output contain instructions for how to update the parameter data and how to read the observation data. The first line uses the PEST template ExtracWells.tpl to write the new AnAqSim file ExtracWells.anaq for the next model run. The second line uses the PEST instruction file ExtracWells.ins to read the model output from run1.out. PEST interfaces with the AnAqSim model using the instructions in these two files.

The DOS batch file used for this example is the command line call of our AnAqSim model.

run_anaq.bat

1	"C:\program files\fitts geosolutions\anaqsim\anaqsim.exe" ExtracWells.anaq run_anaq.txt

This instructs the computer (from the DOS command prompt) to run anaqsim.exe on the file ExtracWells.anaq and execute the commands in run_anaq.txt.

run_anaq.txt

1 2 3

outputfile run1.out calibration exit

The text file run_anaq.txt contains three lines, a line to use run1.out as the outputfile, a line to execute the calibration command so the modeled heads at each of our calibration targets are recorded for each model run, and a line to exit and close the model.

ExtracWells.tpl

The template file is a replica of AnAqSim XML model input file ExtracWells.anaq with the addition of ptf @ at the beginning and each instance of one of our parameters replaced with @<parameter name>   @. The @ symbol (or whatever symbol follows ptf) is used by PEST to identify locations where parameter data should be written. A snippet of the template file is shown below with the relevant sections in blue.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42

ptf @ <?xml version="1.0" standalone="yes"?> <ModelInput> <Release> <ReleaseNo>AnAqSim release 2015-1 beta 30 Jan 2015</ReleaseNo> </Release> <General> <Steady>true</Steady> <Length_unit>feet</Length_unit> <Time_unit>day</Time_unit> <Comments>AnAqSim Remediation Toolkit</Comments> </General> <Solution> <Underrelaxation>0.9</Underrelaxation> <Maximum_iterations>20</Maximum_iterations> <Starting_heads_source>file</Starting_heads_source> <Almost_dry_fraction>0.01</Almost_dry_fraction> <Interface_leakage_option>false</Interface_leakage_option> </Solution> <Check> <Head_check_tolerance>1e-3</Head_check_tolerance> <Q_check_tolerance>1e-3</Q_check_tolerance> <Qn_check_tolerance>1e-4</Qn_check_tolerance> <Extraction_check_tolerance>1e-6</Extraction_check_tolerance> </Check> <ConfinedUnconfinedDomains> <Label_unique>aquifer</Label_unique> <Level>1</Level> <Domain_Type>confined</Domain_Type> <Top_elevation>50</Top_elevation> <Bottom_elevation>0</Bottom_elevation> <Average_head>87</Average_head> <Porosity>0.25</Porosity> <Storativity>0.10</Storativity> <Specific_yield>0.1</Specific_yield> <K1_horizontal>@kh_1 @</K1_horizontal> <K2_horizontal>@kh_1 @</K2_horizontal> <Angle_K1_to_x_axis>0</Angle_K1_to_x_axis> <K3_vertical_top_half>30</K3_vertical_top_half> <K3_vertical_bottom_half>30</K3_vertical_bottom_half> </ConfinedUnconfinedDomains> …

The parameter kh_1 is inserted for both the K1_horizontal and K2_horizontal in the confined aquifer domain for our model. PEST recommends that the variable name be padded by multiple spaces to allow enough space to enter variables (e.g. if the number of desired significant digits for kh_1, above, is less than the space provided, kh_1 will be rounded to fit within the allocated space – 14 digits in the example above).

ExtracWells.ins

1 2 3 4 5 6 7 8

pif @ @Label@ l1 w !cal1! l1 w !cal2! l1 w !cal3! l1 w !cal4! l1 w !cal5! l1 w !cal6!

The instruction file is used by PEST to read any output files and extract the relevant observation data. The @ symbol is used to identify characters in the output file that PEST will place the cursor. The following ExtracWells.ins is used to read run1.out.

ExtracWells.ins begins with pif @ indicating that this is a PEST instruction file and the @ symbol is used to identify characters. The second line, @Label@, lets PEST know that the next command should be run with the cursor placed after the characters Label. In our run1.out file this is the line before the calibration target results. The next line, l1 w !cal1!, indicates that PEST should move down 1 line (l1) place the cursor just after the 1^st white space (w) and read the next uninterrupted set of characters to cal1 (!cal1!). The next 5 lines do the same thing for cal2 through cal6.

run1.out

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

ExtracWells.anaq opened Start solving system 11/20/2018 9:16:42 AM 80 equations and unknowns in system Iteration 1 2 3 4 5 Converged within the specified solution check values. Solve complete 11/20/2018 9:16:42 AM Head Calibration: Label Modeled Head Observed Head Residual Time c-1 89.34256 89.34256 -3.870564E-08 0 c-14 80.57272 80.57272 -4.054981E-06 0 TP-5 97.98874 97.98874 -2.511405E-06 0 MW1 82.09399 82.09399 -1.832374E-06 0 MW-2a 93.70227 93.70227 2.057103E-06 0 TW-84 95.60841 95.60841 2.81206E-06 0 Average residual = -5.94717E-07 Sum of residuals squared = 3.824847E-11 Root mean square error (RMSE) = 2.524826E-06

The final PEST file that may be needed is the parallel run management file, ExtracWells.rmf. This file is used to assign any number of PEST slaves (i.e. the individual processors that will run instances of the AnAqSim model) to run in parallel across a range of locations. For this example problem we used two slaves. These are in the folders slave_1 and slave_2 as subdirectories of the folder where ExtracWells.rmf is located. Each folder that is being used for PEST slaves needs to have a copy of the model, and all of the previous files mentioned. A command prompt then needs to be opened in each of these locations and the command pslave needs to be run to set up these locations to receive instructions from parallel PEST (PPEST) followed by the model command line call (run_anaq.bat in this case). The parallel pest run can then be executed in the main file location by entering ppest ExtracWells into the command line (no extension needs to be provided and the .rmf file and .pst file must share the same name).

ExtracWells.rmf

1 2 3 4 5

prf 2 0 0.2 0 'S1' .\slave_1\ 'S2' .\slave_2\ 1 1

The first line of ExtracWells.rmf indicates that this is a PEST run file, the second line has four entries: the number of slaves (2), whether files are individually named for each slave (1) or not (0), the time to wait for machines to catch up in seconds (0.2) (this should be longer when slaves are distributed over multiple nodes on a crowded network), and whether the lambda search is undertaken using partial parallelization (1) or in serial (0). The third and fourth lines contains the name and location of each slave. For this example the slaves are located in subfolders slave_1 and slave_2. The final line is an estimate of the amount of computation time the model requires to finish on each slave in seconds (if the time is unknown use a ratio of computation times). This is useful when running on nodes of different computational power (e.g. different processors on multiple computers).

PEST Run Results

Pest Run Results

For the box model presented above, PEST performed 68 model runs and was able to complete the optimization with a final phi of 0! The lower the phi value, the better is the estimation of the sought parameters. PEST is telling us that, for this simple ideal model example, the parameter estimation is close to perfect.

The comparison of the calibrated model that was used to set the observations heads and the final PEST calibration is presented below. PEST as able to almost perfectly reproduce the original calibrated model.

Final Calibrated Model

Helpful Tips

PEST utilities are available at http://www.pesthomepage.org/PEST_Utilities.php. These commands can help ensure that the input files are formatted correctly and provide debugging information. The most important one is PESTCHEK <name.pst> which checks that the PEST input files are formatted correctly and consistent.

We hope this simple example has shown you how easy it can be, with a few data and instruction files, to set up and execute a PEST parameter estimation process for an AnAqSim model. If you have questions, please contact us.

The files used for this tutorial are available from FlexAEM at: http://www.flexaem.com/downloads/PEST AnAqSim Example.zip

Exploring the power of AnAqSim using QGIS

FlexAEM recently announced the addition of a new set of exercises to the AnAqSim Instructional Series! The new set of exercises demonstrates how groundwater professionals can use the spatial data management and analysis features of QGIS, a free and powerful open-source geospatial information system (GIS) software package to enhance the groundwater modeling capabilities of AnAqSim. 

The series illustrates how QGIS can generate spatial data describing hydrologic and hydrogeologic features to facilitate groundwater model construction, and then import and process AnAqSim model output data to display model results such as water level contours and flow pathlines in both 2D and 3D. For AnAqSim users who may be less familiar with QGIS, the series includes an overview of QGIS and its capabilities, including compatible file formats, vector and raster analysis tools, editing and management tools, and much more.

The new series includes exercises that demonstrate how to perform a number of useful functions with QGIS, including the digitization of site features to create model elements. For example, river or ocean boundaries, or groundwater divides, can be digitized and transferred to AnAqSim to represent model boundary conditions (Figure 1). The boundaries of other features, such as lakes and ponds, can also be digitized in QGIS and placed in the AnAqSim model using appropriate element properties (Figure 2).

Figure 1

Figure 2

The exercises also demonstrate how QGIS can be used to create rivers and assign elevation data to points along a river system. Users will learn how to digitize a river system using a basemap, and then use a raster file to obtain elevation values for the selected points along the river. Users also learn how to easily format the river data and import it into AnAqSim.

Basemaps for digitizing model elements or preparing river data can also be created in QGIS for use in AnAqSim, or basemap images from other sources can be imported into QGIS to facilitate basemap creation. There is a multitude of freely available GIS data available online, which can be downloaded, viewed, and analyzed in QGIS. The exercises included in the set demonstrate how to use these data to create custom basemaps, and then export those basemaps as a DXF file for use in AnAqSim. For example, the set include an exercise that provides a step-by-step tutorial on how to obtain a web map hosted by the USGS, digitize the boundary of a feature using the web map, and then export the digitized boundary to AnAqSim to help establish model boundary conditions.

Figure 3

Figure 4

Finally, the exercise set demonstrates how to export water level contours (Figure 3) or groundwater flow pathlines (Figure 4) generated from an AnAqSim model run, import the results into QGIS, and view them in 3D. Any model result from AnAqSim that contains elevation values can be viewed in 3D in QGIS, which allows users to easily visualize model scenarios.

To learn more about QGIS and how it can be used with AnAqSim, check out the new addition to the AnAqSim Instructional Series at http://www.flexaem.com/tutorials-tools/Instructional-Set-7.