flexAEM Blog: 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

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.

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

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

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

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

Friday, January 25, 2019

Using PEST with AnAqSim

No comments:

Post a Comment