Viewxsxn reads or generates a set of stacked geologic tops, stored as a set of digital elevation model (DEM) grids, to create 3D solid model reconstructions for use in ground water and other models. In ground water models, the solid models define the 3D character of the ground water aquifers and confining units. Viewxsxn can read and write Surfer, ArcView, or user formatted grids so that geologic tops from other programs can be imported to construct solid models, and imported and/or constructed solid model surfaces can be displayed as structure contour maps. Using Grapher for display, editing, and printing, Viewxsxn can also draw 2D colored cross sections throughout the stack along any user-specified, map orientation in order to visualize the 3D solid model in cross section (xsxn). Viewxsxn can generate a set of top grids based on either projection or folding equations, using a geologic map, a definition of the stratigraphy, and a landsurface DEM for the geologic domain [see From a Geologic Map to a Hydrogeologic Solid Model in Appendix II]. The solid model tops can be incorporated into a Modflow ground water model using additional Geographic Information System (GIS) grids defining river networks / watershed boundaries, aquifer parameters, and recharge. Alternatively, the set of grids comprising the solid model can be output and brought into other model preprocessing software (GMS, Visual Modflow, etc.) to build the solid model for Modflow or other ground water models. A GMS borehole file can be generated from the tops at each grid node or desired multiple of grid nodes. The ability to read and write grids makes the program a useful utility for converting grid formats. To view a cross section, you need:

Golden Software's ( http://www.golden.com/ ) Grapher program.

Installing and Running Viewxsxn

Viewxsxn can be installed manually or using a setup application. Download the setup application, or expand it from a zip file. If you have previous versions of Viewxsxn on your system, you must uninstall the program using the Start\Settings\ControlPanel\Add/Remove Programs routine before installing the new version. Run the setup application by selecting it in Windows Explorer and double clicking the filename. This will install Viewxsxn on your computer, copying 3 program executables to the installation directory and creating an icon on your desktop. Double click the desktop icon to run the program. Should you delete or otherwise lose the desktop icon, launch (or set a shortcut to) the Windows diosxn1.4.exe program (a VB process) in the install directory. It drives the diasxnio.exe and geofold.exe programs (both are Fortran90 processes). The three processes communicate with files, so they are best run in their own directory (which the install routine will setup for you).

User Instructions:

1.0 The Main Menu: The Cross Section Control Panel Form

When Viewxsxn is run, a Multiple Document Interface (MDI) form is displayed on the screen, hosting the Cross Section Control Panel form which in turn hosts the main menu buttons:

At first the cross section starting point, ending point, and zooming (upper z and lower z) textboxes all say "no grid," and all but two of the processing buttons are inactive. At this point, the user chooses between reading or generating a set of geologic top grids for the domain.

1.1 the Open New Set of Grid Files button: reading external grids

When the Open New Set of Grid Files (ONSGF) button is clicked, a message box asks whether the user wants to use an existing “gridfiles.txt” file, located in the Viewxsxn program directory, to load previous sets of grids into the domain. The “gridfiles.txt” file enables you to bypass the input prompts to quickly load pre-existing sets of grids. If you have prepared a “gridfiles.txt” file following its format, or have automatically generated the “gridfiles.txt” file from a previous use of the ONSGF button, and would like to load / re-load these sets of grids, answer “yes.” If not, answer “no”to be prompted for the individual input filenames.

If you have chosen to be prompted for the individual input filenames for each set (answering “no” above), you will be presented with the Sets and Tops form for each set of grids to be read:

The first text field of the form defines the number of sets of grids to be loaded into your domain. YOUR DOMAIN WILL BE DETERMINED BY THE FIRST GRID READ FROM THE FIRST SET, i.e. the first grid will determine the number of rows and columns, the grid spacing, and the coordinates of the lower left and lower right cell. All subsequent grids must correspond to this domain. Reading multiple sets of grids allows you to define different stratigraphy / structures in different sub-regions of your domain. For example, you may use Viewxsxn multiple times to define different structures in different fault-bounded sub-regions, output the gridfiles generated for each sub-region, then re-load all of the sub-regions as patches of the whole domain, each with its own set of gridfiles. The first text field will be grayed out (disabled) for any subsequent sets of grids, displaying the total number of sets to be read as entered on the form for set number 1.

The second text field of the form defines the number of tops (grids) to be opened for the set. For each top, an OpenFile Common Dialogue form will be presented to define the filename(s) of the grid(s) to be imported. The grids must share the same number of columns and rows and the same maximum and minimum x's and y's. The [Surfer or ArcView] header, or the user-specified dimensions, of the first grid of the first set will be used to define the whole domain.

The third text field of the form allows the user to define the ibound integer number for the set of grids. A unique ibound integer thus uniquely defines the corresponding sub-region in the ibound array, generated by Viewxsxn combining subsequent sets into an ibound for the whole domain.

When you click on the OK button, an OpenFile dialogue queries for the filename for each top of each set. The form is a standard windows OpenFile form, but with a check box at the bottom to mark if you want to set, or reset, a blanking value to define regions where the grid is inactive, or to open an ibound (blanking) file for the grid (see below).

If the input file you specify does not have a recognizable extension (*.grd for Surfer or *.asc for ArcView), the Input Format form is presented to specify the “header” values, i.e. the number of columns and rows, the grid spacing, and the minimum x and y values of the lower left node; and the format of the grid:

The minimum (x,y) coordinate specified should correspond the position of the z-value itself (cell centroid for MODFLOW-like, Arc-like, or other block-centered grids; the grid intersection for Surfer-like or other mesh-centered grids). Alternatively to defining the number of columns and rows, the grid spacing, and the minimum x and y values of the lower left node, a convenient sidestep is provided with the IMPORT CURRENT DOMAIN button. If a grid of a domain has already been read, the button will import its “header” values to define the domain. The format is then selected with a radio button, choosing between 1) ascii formatted, where a fortran format is entered in the text box; 2) ascii list directed, i.e. free formatted which is space or comma delimited; 3) fortran's "unformatted," which is a type of binary file containing record information; and 4) binary. The input arrays must be of type real. If using fortran formatted ascii input, be sure to use a real format string. The default format string provided, 6e13.6, is used by Viewxsxn to read and write its Modflow related grids. When the information on the form is complete, select the "read grid with current settings" button.

If you want to 1) reset the Arc or Surfer default blanking value, or 2) set a blanking value for a user-formatted grid, or 3) specify a filename for an ibound array: mark the checkbox, “Set/reset a blanking value or open an ibound (blanking) file for the grid.” CAUTION: blanking values or ibound arrays will blank input at this point. If in doubt about whether you want to control this or not, it is recommended to leave the box unmarked and accept the defaults. Blanking values are used to determine inactive areas, defined as 0 in an integer ibound array that later can be output to a file. Surfer's default blanking value is 1.70141E38, and values greater than or equal to this value are blanked. ArcView's default blanking value is -9999, and values less than or equal to this value are blanked. If none of the values of the grid exceed the blanking limit (Sufer-positively or Arc-negatively), all cells will be active. The ibound array for active cells will be set to the integer specified on the “Sets and Tops” form, and 0 for inactive cells. Alternatively to defining active/inactive cells with blanking values, the ibound array itself can be read in at this point. If you have marked the check box, the Blanking Value or Blanking (Ibound) File form is presented with a radio-button choice between setting a blanking value, or opening an ibound file to be read:

Set the value and select the "set value" button, or select a filename and select open file. If the extension of the file is not recognized (*.grd for Surfer or *.asc for Arc), an Input Format form is presented to specify the format of the integer grid, similar to the Input Format form for reading real arrays described above. An ibound array must be of type integer. If using fortran formatted ascii input, be sure to use an integer format string.

Ibound arrays in MODFLOW are positive integers for grid nodes that are active, 0 for inactive nodes, and negative where values in corresponding arrays are to be fixed. Blanking therefore allows the user of Viewxsxn to make MODFLOW ibound arrays from Surfer or Arc blanked grids.

Once all files are opened, the program will set the domain ranges in the Cross Section Control Panel (main) form.

1.2 the Generate New Set of Grid Files button (see section 2.0)

1.3 the Save Grid File(s) in New Format button

–Initially deactivated. Button activates when tops are read or generated

The stack of grid files can be exported to Surfer, ArcView, or your own specified format files. This is useful for taking Surfer grids (which are "inverted") and writing them for Modflow, or converting between Surfer and Arc formats. The conversion between Surfer and Arc grids will also convert the blanking values between each. Surfer grids report a min and max z-value for the grid in the header. When outputting the stack to individual Surfer grids, the maximum z and minimum z for each file is recalculated (1 set for each layer). These are different than the minimum and maximum z's that are reported for the whole stack on the main form.

When the button is selected, a Save File Common Dialogue (Save Top As) form is presented; the form additionally has a checkbox option for setting output blanking values with a textbox for the value, and a checkbox option for outputting the existing ibound array, determined from either reading or generating the grid:

If the “Reset Blanking Value” checkbox option is selected, enter an output blanking value in the textbox. The blanking values entered here do NOT determine inactive regions in the output at this point, but are written into the output where inactive areas have already been determined. Inactive areas are originally determined 1) on input (by input blanking values), or 2) an input ibound array, or 3) by ibound arrays produced in generating grids (see 2.0).

1.4 the Read Modflow96 heads from Binary File button

–Initially deactivated. Button activates when tops are read or generated

If a Modflow96 model has been run and a binary head output file produced, this button will read the binary heads and convert them to ascii grids for each layer for each time-step for each period. To make a binary head output file, use the the Make Modflow96 Model from Tops button (see section 1.5 and 3.0). The Read Modflow96 Heads button is not available until tops have been read or generated. If running Viewxsxn fresh, and/or the domain of the Modflow model has not been entered into Viewxsxn, it is suggested to read the DEM of the model domain to activate the button as well as define the domain. When the button is selected, the user is prompted to provide the Modflow96-defined NAM file, which lists the input and output files associated with a Modflow96 run. By reading the NAM file, Viewxsxn determines the simulation’s dimensions and time-steps, and reads the NAM-file-specified binary output file. It then writes one grid for each layer for each time-step for each period. The output filenames are:

head#_#T_#L.txt

where # is the period number, #T is the time-step within the period, and #L is the layer number. The output format is 6e13.6.

These head-text files can be read using Viewxsxn’s Open New Set of Grid Files button (see section 1.1) and then converted to Surfer, Arc or other formats using the Save Grid File(s) in New Format button (see section 1.3). Because of the “txt” file extension, reading these files will result in the form for domain and input format specification (see section 1.1). If re-running Viewxsxn and/or the domain for the Modflow files has not yet been defined, one suggestion is to read the corresponding DEM domain first, then select the “Open files” button again and read the head-text files, using the IMPORT CURRENT DOMAIN button to import the DEM definition. Then select the format radio button: if the head-text files were generated by Viewxsxn, the format of the grids is already specified, 6e13.6.

1.5 the Make Modflow96 Model from Tops button (See section 3)

–Initially deactivated. Button activates when tops are read or generated

1.6 the Write GMS Borehole File button

–Initially deactivated. Button activates when tops are read or generated

The Write GMS Borehole File button enables the user to output a borehole file (of "artificial boreholes" ) formatted for input into GMS. Grid information constructed for domains of DEM's is commonly too dense to be used in the solid body module of GMS. The ability to write a borehole file enables the Viewxsxn user to take a smaller sampling of the data for input into GMS. This routine (button) writes a lattice of boreholes, up to the number of columns in the row direction (column to column) and up to the number of rows in the column direction (row to row). For each borehole location, all the tops from the grid stack that exist at that point (i.e. have non-zero ibounds defined) are written to a borehole. When the button is selected, the following saveas command dialogue form is presented:

The form has two additional textboxes, one for the number of boreholes in each (row and column) direction. The program writes the number of boreholes in each direction equally spaced, and only writes a borehole if it is within the boundary of the problem (where ibound is not zero). Therefore it is possible to select for, say, 30 boreholes along a row, but have varying numbers of boreholes less than or equal to 30 in each row depending on the shape of the domain boundary. The program "boundary checks." It will not let you specify a number of boreholes in the row direction greater than the number of columns, or in a column direction greater than the number of rows. If the numbers are greater than the number of columns or rows respectively, the program sets the value to the number of columns or rows respectively. If the numbers are less than 1, the program sets the values to 1. Select the filename for the output and select the Save Boreholes to File button.

Note: Though boreholes efficiently write the tops for a single location, it is not efficient to write a large number of boreholes; a large number of boreholes results in a very large output file. For instance if the number in each direction is set to the dimensions of the grid, the file will be much larger (in bytes) than the sum of the sizes (in bytes) of each grid file comprising the tops. This is due to the extensive header information that must be written for each borehole. For example, a test domain that was 361 cols x 470 rows x 8 layers required about 1.5M for each of 8 grids (for a total of about 12M) whereas the borehole file was 54M.

1.7 the Set/Reset Colors and View buttons (Viewing a Cross Section)

–Initially deactivated. Button activates when tops are read or generated

To view a cross section, 1) set the coordinates of the starting (left side) point and ending (right side) point of the cross section, 2) set the upper and lower z value of the view (for zooming), 3) select the "set/reset colors" button if you want to set or reset the colors of the cross section, and 4) select the "View" button. Note: selecting the View button requires that Grapher ( www.golden.com ) be installed to view the cross section automatically. Cross section data files are generated that can be viewed with other programs.

Coordinates can be specified in either the column-row coordinate (i,j) or the x-y coordinate. The program "boundary checks." It will not let you specify an i,j or x,y that is out of range. If an x or y, or an i or j, is too low, it will default to the minimum value allowed. If it is too high, it will default to the maximum value allowed. If you specify a column (i) coordinate, the program will calculate the x coordinate for that column and will update the associated text field. Same for j and y. If you specify an x coordinate, the program will calculate the i and will update the associated text field. Same for y and j. The program can thus be used as a grid calculator to go between i,j and x,y coordinates.

The upper and lower z value will determine the upper and lower limits of the cross-section in grapher. It is defaulted to the minimum and maximum z determined for the whole set of grids. The minimum and maximum z defaults can be restored with the "min" and "max" buttons.

Cross sections initially do not have color unless the "set/reset colors" button has been selected and the routine has been completed. A color palette for setting the layer colors is presented; one palette is presented for each layer in the stack. For the grid nodes used in a colored cross section, the color specified for the layer is used between the layer's z-value (considered a top) and the next layer's z-value at any given node. Once set, colors will remain unchanged unless the user changes the number of layers either read or generated, or resets the colors, and selects the View button again . If the number of layers is changed, the default colors are set to white (no color).

1.8 the About button

The current version and other program information is displayed with the About button:

1.9 the Quit button

The Quit button terminates the program and its processes.

2.0 the Generate New Set of Grid Files button: producing a section using projection or folding

When the Generate New Set of Grid Files button is selected, the Fold and Projection Parameters form is presented, tiled with the main form in the multiple document interface:

All options on the Fold and Projection Parameters form are initially inactive, except for the Open New Geologic Map and Elevation Files button.

2.1 the Open New Geologic Map and Elevation Files button

Selecting the Open New Geologic Map and Elevation Files button will present an OpenFile Common Dialogue form for opening a geologic map, followed by another to open the elevation file (DEM) of the same domain. Both are required for Viewxsxn to generate a stack of geologic top grids to form the solid model. Again, Surfer, ArcView, and user-formatted grids are supported on input using the same procedure described in section 1.1. When opening the geologic map file, the standard Surfer or Arc files blanking values (from either the Arc header or the Surfer default blanking value), will be used to define inactive areas of the domain. The checkbox at the bottom of the OpenFile dialogue form is used to either set/reset the blanking value, or import an ibound array for determining inactive areas of the map, as described in section 1.1. The header-defined blanking value, or user-defined blanking value, or imported ibound array define the inactive (and thus active) portions of the domain for the calculation of a set of projected or folded horizons.

The geologic map file must be an integer grid (array) with unique integers for each mapped horizon to be projected or folded. Think of it as geology-by-numbers. The map can easily be generated in ArcView by gridding geologic map shapefiles (see section Appendix II below). The grid domain must correspond to the same domain as the elevation file. The elevation file can contain real or integer values (they will be processed as real numbers).

When the reading of the geologic map and elevation file are complete, the Stratigraphic Units and Top Assignments (SUTA) textbox and the Number of Tops to Generate textbox are activated, and the integers encountered in the geologic map (geology-by-numbers) are listed in the SUTA textbox.

2.2 the Number of Tops to Generate textbox

In the Number of Tops to Generate textbox, enter the number of tops (grids) to be generated. This number (nz) can be less than, equal to, or greater than the total number of geologic units (nu) found in the geologic map (geology-by-numbers) and listed in the SUTA textbox. Any extra tops are appended to the bottom of lowest mapped unit (see sxn 2.4 below).

2.3 the Stratigraphic Units and Top Assignments

After the geologic map and elevations have been successfully read, the integers encountered in the geology-by-numbers map are listed in the Stratigraphic Units and Top Assignments (SUTA) textbox. These integers may or may not be numbered in stratigraphic order; they are just the stratigraphic units (in integers) of the map. They need to be given top assignments in stratigraphic order. A top assignment must be provided for each stratigraphic unit by clicking behind the number on each line, entering at least one space, and entering a top number. The top numbers must be ordered stratigraphically downward, with 1 being the topmost top and nz (the total number of tops) being the lowermost layer. For an example, see Appendix II. Numbering downward is consistent with the MODFLOW convention for numbering layers. The same top number can be assigned to multiple stratigraphic unit numbers listed in the SUTA textbox; this results in stratigraphically lumping units on the map to one unit with one top assignment.

If the geology-by-numbers map was originally constructed (e.g. gridded from ArcView) with integers that are already downward-stratigraphically-ordered as desired, then all one needs to do is repeat each number, i.e. assign the same top number as the stratigraphic unit number, on each line. Note about ArcView: unless the gridding was based on an explicitly defined field of stratigraphically ordered numbers for each record, the gridding will produce random numbers that cannot be assumed to be in stratigraphic order, even though the numbers may be sequential.

Lumping: If it is desired to lump two or more successive stratigraphic units into a single unit, assign the same top number to each of the stratigraphic units to be lumped. If the total number of tops (nz) is less than the total number of geologic units (nu), i.e. nz < nu, the extra stratigraphic units will be assigned to the top that defines the top of the group of lumped strata. If lumping, be sure to set the number-of-tops-to-generate (nz) equal to the total number of tops assigned, which will be less than the total number of stratigraphic units listed in the SUTA. If the number-of-tops-to-generate exceeds the number of top assignments made, extra tops will be appended to the bottom of your map-defined stack (see next paragraph on “Appending”).

Appending: If the total number of tops is greater than the total number of geologic units (nz > nu), or if only a subset of total number of tops (nz) is actually used in making top assignments (the remaining tops combined into a top assignment), the additional tops will be appended to the bottom of the stratigraphically-lowest top assigned from the map.

2.4 the Strike and Dip text boxes

In the strike and dip textboxes, enter the average strike and dip that will be used for projection, or is most representative of the outcrop widths on the geologic map in the region of the folding. For folding, these outcrop widths and the strike and dip will be used to calculate "depth thicknesses" (the vertical thickness of the strata as opposed to the actual stratigraphic thickness). See the discussion in Appendix I for details on projection and/or fold thickness calculations.

2.5 the Calculate Vertical Depths button and Vertical Depths to Next Top textbox.

After setting the number of tops to generate (section 2.2), making the stratigraphic assignments (section 2.3), and providing a strike and dip (section 2.4), click the Calculate Vertical Depths button to determine the vertical “thickness” between the tops. A thickness from the top to the next top downward will be calculated and listed for each top in the Vertical Depths to Next Top textbox. The thicknesses are based on contact widths from the “geology-by-numbers” map, and are a function of the combination of the number of tops, stratigraphic assignments, and dip provided (see Appendix I). An average for the whole domain (map area) is reported. At this point, the user can change the number of tops, stratigraphic assignments, and dip experimentally until attaining reasonable vertical “thicknesses” between units. The first top, and any other top that does not have BOTH an upper and lower contact defined on the map will report a zero thickness. If appending unmapped tops to the bottom of the lowermost mapped unit (see section 2.3), zero thicknesses will be reported for each. The last top will list zero thickness since there is no top defined below it.

Once a set of thicknesses has been calculated, the user can either accept the given thicknesses, or override by typing new values in the text field. The thicknesses are only used in folding reconstructions (see below), and for appending units to the bottom of the lowermost mapped unit in both folding and projection reconstructions. The final values in the textbox when initiating these reconstructions will be the ones that are used. Though the textbox thicknesses are only used for optionally appending the bottom of projection reconstruction, the thicknesses calculated and listed in the text give a sense of the thicknesses that will result in projection using the dip value.

2.6 the projection option (radio button)

If this option is selected, each layer will be projected from its contact using the formulas in Appendix I. Viewxsxn assumes that the minimum distance to the map contact of the horizon is along a traverse perpendicular to strike ( i.e. parallel to the dip direction). Because the minimum distance to a curve is perpendicular to the curve at the point of intersection, this essentially assumes that, on average, the map contacts run parallel to strike. In most cases, particularly in folded regions with high dip angles, this assumption is justified. However in some cases, a contact may run parallel to dip, particularly along slopes that are unrelated to the structural trend, and this assumption is not valid. After several experiments with different actual geologic maps, this assumption resulted in better projections than projections corrected for apparent dips from an assumed fixed strike direction, primarily because the strike (and thus dip direction) changes along most structures.

Though the projection option is mathematically more simple (simple planar projection) than the folding option, it is computationally more intensive because the formulas depend upon the calculation of the minimum distance between the point and a point on the contact. The contact is stored in the raster of the whole domain (rather than in a vector curve), and thus the entire raster must be searched for each point in the subcrop/outcrop region where the grid is calculated. When the "Calculate Tops Using Current Settings" button is selected, the program may appear to "lock-up." If the lock-up is truly not computer related (or a program bug), wait long enough for the calculation to finish, and control of the program will be restored.

2.7 the folding option (radio button)

If this option is selected, the folding-parameter textboxes will become activated. Fill out the parameters for the fold and select the “Calculate Layers Using Current Settings” button.

A fold horizon is generated for the lower-most (highest numbered) layer defined. Successive (lower numbered) layers are calculated by adding the “depth thicknesses” (see Appendix I) to the basal fold horizon. Both folding and plunging are periodic in the Valley and Ridge. Though double plunging is most likely controlled by relief along buried thrust ramp surfaces, a regular pattern of plunging suggests periodicity. The fold horizon is calculated by orthogonal sine waves, one describing folding (in the y-direction with axial traces in the x-direction) and the other describing plunging (in the x-direction):

(1)

with wave numbers

(2)

and

(3)

where λ₁ and λ₂ are the wavelength of the folding and the plunging respectively, A and B are the amplitude of the folding and the plunging respectively, and x₀, y₀, and z₀ are offsets for the origin. The user inputs A, B, x₀, y₀, z₀, and wavelengths λ₁ (l1 on the form) and λ₂ (l2 on the form) in the textboxes on the form. In addition to the offsets for the origin providing arbitrary translation, an angle of rotation is calculated from the strike to provide arbitrary rotation for the coordinate system. The result is that the axial trace of the fold (the x-axis of the rotated coordinate system) will parallel the strike provided in the strike textbox.

Note: Strike is assumed to be measured clockwise from the positive y-axis of the grid which is assumed to parallel north. If your grid is in some other orientation, just add or subtract the difference from the strike direction. For instance, if the grid’s y axis trends N30°W and the strike is N45°E, the strike provided in the textbox should be 75°.

Section 3.0 the Make Modflow96 Model from Tops button: Generating a Modflow Model.

To demonstrate, the Mahantango field example, featured in Appendix 2 to demonstrate solid model construction, is further used here to demonstrate the construction of a Modflow model from the geologic tops in Viewxsxn. Once geologic tops have either been generated or imported into Viewxsxn, the resulting solid model can be incorporated into a Modflow96 model. The Make Modlfow96 Model from Tops button initiates a Modflow collator, i.e. a routine that will collate the constructed solid model and other previously-defined Modflow parameters into Modflow96 input files. Although Viewxsxn can construct an entire Modflow96 model on its own, the collating routine is not a pre-processor in the strictest sense. Compared to commercially available pre-processors, it provides limited tools for model creation, and currently does not accommodate Modflow model / input-file editing. Thus the collating routine requires that the parameter definitions be ready before-hand, whether they are single values or arrays, and information generated from other pre-processors can be helpful in this regard. If using Modflow’s River Package, setup in Viewxsxn requires an integer-raster, river-definition file, coincident with the solid model domain (see section 3.x). This file can usually be created in ArcView (see Appendix II on raster output of line and shapefiles).

When defining parameters to be used in Modflow, dimensional units are your business. Modflow requires all units of length to be consistent, including those in definitions of pump rates, conductivities, cell dimensions, surface elevations, transmissivities, and specific storage coefficients, etc. Modflow also requires the user to pick a time unit, which will be used or assumed in all rate terms. Viewxsxn does not convert units. In short, the routine launched from the Make Modflow96 Model from Tops button is a Modflow collator for the Modflow initiated.

3.1 the Modflow Setup form

When the Make Modflow96 Model from Tops button is selected, the Modflow Setup form is displayed.

At this point, the user defines the number of layers in the Modflow model to be defined from the tops in Viewxsxn. Next the user selects whether the Modflow simulation to be defined is transient or steady state, and the time unit that will be used in all rate terms, periods, and time steps. Next the user selects the Modflow packages to use. Viewxsxn automatically builds the Basic, Block Centered Flow, and Output Control packages. With this version of Viewxsxn, the user can optionally add the River and Recharge packages. For the Mahantango example, a 4-layer, steady state model using seconds for time, and employing the River and Recharge package is shown in the form.

In the last field, the user chooses a method to handle pinchouts, and whether or not to equate parameters in model layers that are defining the same aquifer in pinchout zones. “Highly recommended” options are preset as defaults. When constructing a model, the user is encouraged to thoroughly read section 3.2 before changing these options. Finally, a textbox is provided to set a minimum thickness for layers. This value will define minimum layer thicknesses in both pinchout zones and in regions where imported model surfaces may be ill-defined (bottoms defined above tops, etc.). An extensive discussion of these pinchout options using the Mahantango example is discussed in section 3.2.

When the form is filled out, selecting the Setup Modflow Model button will launch the remainder of the model setup, starting with a Laycon and Top/Bot Assignments for Modflow Layer# form for each model layer defined.

3.2 the Laycon and Top/Bot Assignments for Modflow Layer# form

The Laycon and Top/Bot Assignments for Modflow Layer# form allows the user to assign the top surfaces, defined in Viewxsxn that were either generated from the geologic map (Section 2.0) or imported (Section 1.1), to Modflow layer top and bottom surfaces. Then the user defines whether the Modflow aquifer layer is confined, unconfined, or convertible by setting the Laycon value for the layer (please see the Modflow documentation for information on Laycon):

Modflow defines aquifer and confining units using layers, and though these are related to geologic tops, they are not the same thing. Basically, a layer is defined with both a top and a bottom, whereas a top is just a single planar horizon in the subsurface. A layer’s top and bottom are defined with tops. Another key distinction is that a top is just a surface marking the occurrence of a specific geologic unit (rock, stratum, etc.) below it, whereas a layer attempts to define the extent and geometry of a user-defined unit. A top is usually thought of as a geologic contact between units in the subsurface, however it can also be an unconformity surface and a fault surface. Depending on the definition of a unit ascribed to a layer, a layer may have other non-defining tops within, in addition to the tops defining its top and bottom boundaries. With this form, Viewxsxn allows the user to define each Modflow layer’s top and bottom, using the tops reconstructed in, or imported into, Viewxsxn. The user assigns a top number, i.e. a number from the set of tops numbered from 1 at the surface downward to the total number of tops, constructed with, or imported into Viewxsxn, to each numbered layer’s top and bottom, i.e. the Modflow layer numbers also numbered from top to bottom. To accomodate any extra tops within a layer, Viewxsxn allows a user to lump tops within a layer simply by skipping any lumped tops in layer top and bottom assignments.

For the 4-layer Mahantango example, the first Modflow layer number is assigned Viewxsxn top numbers 3 and 4 for the top and bottom respectively. The remaining layers were assigned tops and bottoms respectively as follows: tops 6 and 8 for layer 2, tops 8 and 9 for layer 3, and tops 9 and 10 for layer 4. Please see the pinchout discussion below for the implications of skipped layers and pinchouts with these assignments.

3.3 Pinchout handling:

Most commonly, a top surface of a specific geologic unit will transition in the subsurface from a contact surface between two geologic units into an erosional surface called an unconformity. If the erosional surface cuts across multiple geologic units, the erosional surface becomes a collection of many top surfaces as each contact top surface emerges from below and merges into the erosional top surface at an angle. The cross cutting erosional surface is known as an angular unconformity because of the angle formed between the contact surfaces and the erosional surface. As multiple top surfaces, including layer-defining surfaces, merge into the angular unconformity from below, the tops and bottoms of layers merge together and pinch out along the angular unconformity surface. In the strictest geological definition, the term angular unconformity is reserved for an erosional surface where the contacts above and below the surface are at an angle to each other due to tectonic displacements, i.e folding or faulting, that affected the units below the surface prior to the genesis of the units above the surface. More generally defined as above, the angular unconformity is by far the most common cause of layer pinch-out.

Layer pinchouts are problematic in Modflow and most other ground water models. Most models require layers to be continuously defined throughout the domain. A further requirement of most ground water models is for the mesh of cells defining layers to be “logically rectangular [or logically triangular]; it must be possible (hypothetically) to reposition the nodes to form a regular mesh of cube- [or triangularly prismatic-] shaped elements without adding or deleting any connections between nodes.” (Voss, p. 55, insertions are this author’s extension of the definition). Thus a layer of cells must extend beyond any pinchout. The problem becomes what to do with the layer beyond the pinchout. It is possible to de-activate the layer in Modflow in regions beyond the pinchout by setting the ibound to zero for that layer. However this usually results in layers being defined below layers that are turned off. This is not allowed in most ground water models including Modflow.

Two approaches to handling this model grid-design problem in Modflow have been summarized by Jones et al. (2002), the grid overlay method and the boundary-matching approach. The grid overlay method imposes a uniform grid that does not necessarily match hydrostratigraphic boundaries, but are orthogonal in both the horizontal and vertical plane. The method usually results in strata-cross-cutting “layers” in the model, thus the layers usually define different hydrostratigraphic units, or even combinations of units, in different zones of the domain. Alternatively, the boundary matching approach “ensures[s] that each upper and lower boundary [of the hydrostratigraphy] is precisely matched by a layer boundary in the Modflow grid.” Though perhaps a matter of modeling opinion, there are distinct advantages in using the boundary-matching approach, namely that layers correspond to aquifer units, allowing for a uniform definition of aquifer properties and head reporting within a model layer, and that anisotropy effects can be imposed numerically along unit boundaries (Hoaglund and Pollard, 2003). However, the main disadvantage to using the boundary matching approach is that the imposed coincidence of layer tops and bottoms with geologic tops result in pinchouts along angular unconformities.

The angular unconformity pinchout problem is a major obstacle in defining solid models for regional flow models using geological maps. As stratigraphic entities, unconformity surfaces, including angular unconformity surfaces, are usually thought of as buried erosional surfaces. However, the most common and pervasive angular unconformity is the land surface. Thus a geologic map, loosely defined as a map of the outcrop extents of a given set of geologic tops for a given domain, can alternatively be defined as a map of layer pinchouts at the surface for a given set of layers bounded by a given set of geologic tops for a given domain.

Denuded synclinal and basinal structures and their analogs occasion the most common, problematic model geometry when constructing solid models across the landsurface angular unconformity using geological maps. Moving outward along the unconformity away from the synclinal axis or basinal center, successively older units, deeper within the structure’s core, pinch out leaving successively wider layer extents below layers with lesser extent. Most ground water models including Modflow do not allow active layers below inactive layers. The grid overlay method and the boundary matching approach each produce its own set of problems in handling model layer definition. If employing the grid overlay method, each model layer must have its own set of different zones defining the different extents of each aquifer within the layer. If employing the boundary matching approach, layer pinchouts leave model nodes beyond the pinchout that cannot be deactivated above active layer nodes.

In most ground water models, model layers with wider extents cannot be below model layers with lesser extents; the inverse is not true. A modeling metaphor to help recall this modeling rule: if model layers are layers in a wedding cake, up-side-down wedding cakes are o.k., but right-side-up wedding cakes are not allowed. Denuded synclines or basins along the landsurface angular unconformity are problematic, right-side-up wedding cakes. Ironically though, the problem area associated with the problematic synclinal or basinal structure is in the region of the anticlinal or domal core, where less extensive upper layers do not exist, but where nodes previously used to define a now-pinched-out upper layers cannot be deactivated: the upside-down wedding cake doesn’t exist in the anticline, but rather an upside-down wedding cake set of nodes must be completed to maintain layer continuity.

In order to employ the boundary matching approach with its distinct advantages discussed above, the required redefinition of pinched-out nodes can be handled in Viewxsxn using one of two methods: draping the nodes over the outcrops of the lower layers, or distributing the nodes within the outcrops of the lower layers. “Draping the nodes” means to define layers of unit thickness over the outcrop and then equating aquifer parameters to those of the outcropping layer. A similar method was employed in Hoaglund (2002) where thin layers were draped over subcrops of basinal bedrock below an angular unconformity overlain by glacial drift. A slight difference was that the thin layers were given the larger hydraulic parameters of the overlying drift assuming the layers themselves were negligible. This method is not recommended; though the layers are negligible, it is not clear what the effect of the layers actually is in model calibration. The highly recommended procedure for handling the pinchout problem with Viewxsxn is the option to distribute extra, pinched-out nodes into the outcropping layer.

Figure 1 shows a block diagram generated with Viewxsxn from the geological map, the surface view, and projection reconstruction drawn with Viewxsxn using Grapher, the side view.

The reconstruction corresponds to the Mahantango 4-layer Modflow model example in Appendix II (see section II4), with the cutaway along column 139. and is used to demonstrate the method of distributing Modflow nodes within the outcrops of the lower layers. Geologic top surfaces generated in Viewxsxn are shown numbered 1 to 10 in square boxes drawn on the top. Note that the tenth top does not outcrop, and instead was appended to the bottom of top number 9, and that top 1 has a lower contact with top 2 but not an upper contact. Modlfow layer numbers are shown numbered 1 to 4, with layer-centered nodes numbered 1 to 4 in circles.

Geologic tops were assigned to the 4 Modflow layer’s top and bottom definitions. Geologic tops 1 and 2 are the tops of the Llewellyn Formation (Pl) and the Pottsville Conglomerate (Pp) respectively, a massive section of ridge forming siltstones, sandstones, conglomerates, and coal. These units were hydrostratigraphically interpreted as confining units. Geologic top 3 is the top of the Mauch Chunk Formation (Mmc), the top of the first valley-floor aquifer represented by Modflow layer 1. Geologic tops 4 and 5 are the tops of the Pocono Sandstone (Mp) and Spechty Kopf Formation MDsk, ridge-forming siltstones, sandstones, and conglomerates; geologic top 4 was assigned to the bottom of Modflow layer 1. Geologic tops 6 and 7 are the tops of the Duncannon Formation (MDcd) and Sherman Creek Member of the Catskill (Dcsc) respectively, two geologic formations hydrostratigraphically lumped into one aquifer unit represented by Modflow layer number 2. Geologic top 6 was assigned to the top of Modflow layer 2. Geologic tops 8 and 9 are the tops of the Irish Valley Member of the Catskill (Dciv) and the Trimmers Rock Sandstone (Dtr), forming the top and bottom respectively of an aquifer unit, Modflow layer 3. Geologic top 8 was shared: it was assigned to both the bottom of Modflow layer 2 and the top of Modflow layer 3. Geologic top 9, along with top 10, the top of the Mahantango Formation (Dm), form the top and bottom respectively of an aquifer unit, Modflow layer 4. Geologic top 9 was shared: it was assigned to both the bottom of Modflow layer 3 and the top of Modflow layer 4.

Tracking the Modflow layers from the synclinal to the anticlinal axis right to left along the cross section, the pinchout problem becomes immediately apparent. The mapped positions of the Modflow nodes are shown as a projected grid on the map surface. The actual positions of the nodes in Modflow are layer-centered at depth, as shown along the profile. In the right-most (southern-most) column, layer 1's node is shown centered between its top and bottom, tops 3 and 4. The pattern is repeated at the second next node to the left (north), but by the third node, the Pottsville has pinched out as its bottom, geologic top 3 reaches the landsurface, and the top of the Modflow aquifer layer 1 becomes the landsurface. At this point, in the outcrop area of Modflow layer 1, the layer 1 centered node is half-way between the landsurface and geologic top 4, the bottom of layer 1. As geologic top 4 approaches the landsurface, the node is centered in progressively less thick aquifer unit 1, until the aquifer unit completely pinches out and the layer 1 node must be re-assigned. At that point, the model layer is re-defined, joining layer 2, and nodes 1 and 2 then evenly divide aquifer unit 2 between tops 6 and 8. Note that before the pinchout, layer 2's node solely defined aquifer unit 2, midway between tops 6 and 8. As geologic top 8 approaches the landsurface, nodes 1 and 2 evenly define unit 2, layer 3 defines aquifer unit 3, and layer 4 defines aquifer unit 4. When aquifer unit 2 pinches out, nodes 1,2, and 3 evenly divide aquifer unit 3, and node 4 soley defines aquifer unit 4. Finally, all 4 nodes evenly divide aquifer unit 4 where it is the sole aquifer unit in the core of the anticline. One disadvantage of re-distributing the nodes evenly in the outcropping aquifer unit is that as the aquifer unit thins, the model layers thin and edge effects can be generated in the model head solution along aquifer contacts. This can be moderated by setting the minimum layer thickness correction high enough. A distinct advantage of re-distributing the nodes in the shallowest unit is that this redefinition provides extra nodes in the shallower, more dynamic part of the ground water system.

3.4 The Nam file

After Modflow top and bottom assignments have been made, Viewxsxn ques the user to provide a NAM filename. The NAM file lists and categorizes all of the Modflow model input files into the run-directory. The path provided in the NAM filename becomes the run-directory, and the path to all of the input files and eventual output files.

3.5 Aquifer starting heads, parameter definitions, and recharge values.

Viewxsxn next ques for all of the aquifer parameters using a form with a radio button that allows the user to choose between defining the whole domain with a single value, or providing a filename to an array with spatially varying values. An example is given for the queuing of the vertical conductivity of an intermediary layer:

In addition to the horizontal hydraulic conductivity defined in Modflow, Viewxsxn ques for vertical hydraulic conductivities for all of the aquifer units and any intermediary confining units between layers. Viewxsxn calculates Modflow’s VCONT using equation 52 of the Modflow documentation (McDonald and Harbaugh, 1988, p. 5-16 ) with the definitions of layer tops and bottoms, and vertical conductivities.

After the parameters are specified, Viewxsxn builds the BCF file and writes the model input files to the run-directory, listing each filename in the NAM file.

3.6 The River Package

If the River Package option has been chosen with the checkbox on the Modflow Setup Form (see section 3.1), the user will be prompted for the river raster file. The river raster is a grid corresponding to the model domain where non-zero integers are used to define river cells and zero’s denote interfluves (non-river cells). The program then reads the entire river raster grid, sets the stage of the river cell to the DEM value, and identifies all of the unique integers in the grid, presenting a list of those integers on the River Parameters form:

For each unique integer in the river raster, a unique set of parameters may be defined. These are the layer hosting the river cell, the river depth, the thickness of bottom sediments, the hydraulic conductivity of bottom sediments, and the area of the river reach. The total of the depth of the river plus thickness of the river sediments is subtracted from the stage to determine RBOT. The grid spacing from the domain, provided on the form, can be used as a proxy for the length of the river reach in the cell. The area can thus be determined by multipling the grid spacing times the river width, i.e. the area of the bottom of the river channel, LxW. The area is then used along with the thickness of bottom sediments (M) and their hydraulic conductivities (K) to calculate the Modflow River Conductance at for each set of river integers, using the Modflow formula:

Cond = KxLxW / M

The following figure from the Modflow documentation (McDonald and Harbaugh, 1988) will help the user conceptualize these parameters, but the user is encouraged to read the documentation.

Values are entered by clicking behind the integer in each field, typing at least one space, and then entering the value representative for that integer. The unique integers provide a way to vary river parameters by proxy, e.g. by stream order, or geology, etc. When the form is filled out, select the Write Modflow River File button.

If the recharge package has been selected, the program asks for the recharge value or a file of a grid of recharge values, and the Modflow setup routine is complete.

Modflow can be run by placing a Modflow executable in the run directory, typing the program name on a command line, and answering Modflow’s query for the NAM filename. If Modflow successfully runs and produces a head solution, a binary output file will be generated. These heads can be post-processed into Surfer, ArcView, or other using Viewxsxn (see Section 1.4).

References

Billings, M.P., 1972. Structural Geology 3rd Edition. Prentice Hall. Englewood Cliffs, New Jersey. 606 p.

Burton, W.C., L.N. Plummer, E. Busenberg, B.D. Lindsey, and W.J. Gburek. 2002. Influence of fracture anisotropy on ground water ages and chemistry, Valley and Ridge Province, Pennsylvania. Ground Water 40 no. 3: 242-257.

Hoaglund, J.R., G.C. Huffmann, and N.G. Grannemann. 2002. Michigan Basin regional ground-water flow: discharge to three Great Lakes. Ground Water, v.40, no. 4: p.390-405.

Jones, N.L., T.J. Budge, A.M. Lemon, and A.K. Zundel. 2002. Generating MODFLOW grids from boundary representation solid models. Ground Water 40. No. 2: 194-200.

McDonald, M.G., and A.W. Harbaugh, 1988. A modular three-dimensional finite-difference ground-water flow model. USGS Techniques of Water Resources Investigations, Book 6, Chapter A1, Modeling Technologies.

Stabler, C. L. 1968. Simplified Fourier analysis of fold shapes. Tectonophysics 6: 343-350.

Tarboton, D. G., (1997), "A New Method for the Determination of Flow Directions and Contributing Areas in Grid Digital Elevation Models," Water Resources Research, 33(2): 309-319.

and website:

http://hydrology.neng.usu.edu/taudem/

Voss, C. I., and A.M. Provost. 2002. SUTRA, A model for saturated-unsaturated variable-density ground-water flow with solute or energy transport. U.S. Geological Survey Water-Resources Investigations Report 02-4231, 250 p.

Appendix I. Calculating the depth to the horizons using projection, or calculating the “depth thickness” for folding.

From the given information the program calculates the depth to the horizon during projection, or calculates a “depth thickness” for each horizon to be stacked on the folded horizon. The given information includes: 1) the Digital Elevation Model (DEM): Terrain elevation, z, known every 90 meters, in x and y, or better, 2) Digital Geologic Map: ArcView shape files of local geology converted to an integer raster (the geology-by-numbers grid) at the DEM resolution, 3) (Hydro)Stratigraphic relationships, and 4) the strike and dip measured at outcrops or determined from boreholes. For calculating the depth to a horizon and/or the “depth thickness,” the following formulas are used (for the case where the beds dip in the downslope direction):

If the depth to the horizon is calculated at a point across the outcrop at the next stratigraphically higher contact, the depth becomes the “depth thickness.” The stratigraphic thickness (t) can be calculated from the “depth thickness” (d) by a simple trigonometric relationship [ t = d cos( δ ) ], or it can be calculated directly from the given information using the following formula (for the case where the beds dip in the downslope direction):

If the beds dip in the upslope direction, the formulas become:

Appendix II. From a Geologic Map to a Hydrogeologic Solid Model: An example using ArcView, Viewxsxn (with Grapher), and GMS

Regional ground water modelling requires the definition of aquifer geometries beyond the typical extent of consulting and/or water supply projects where detailed, professional geologic logs may exist. The only data that exists over broad regions are water well records. The lithologic logs of such records are usually [pathetic] insufficient to characterize the geology over the desired regional extent. As an example, ten to twenty separate formations have been identified in the Nittany Valley, a typical carbonate valley in the Valley and Ridge, but most water well records simply define “limestone” from the top of bedrock to the bottom of the well. Yet the differences in physical properties between the different “limestone” formations–particularly differences in the ratio of limestone to dolomite which affects dissolution–causes the contrasts in hydraulic properties, and thus the definition and characterization of aquifer and confining units. Fortunately excellent geologic maps are available for most places in Pennsylvania, and the information from these maps can be used to reconstruct the geology at least to a shallow, but sufficient enough extent for most ground water analyses. The following example outlines a method to generate a solid model (the aquifer geometry for a hydrogeologic model) from a geologic map using ArcView, Viewxsxn (with Grapher), and GMS. To follow this example, the following extensions, tools, and scripts are needed in ArcView 3.2:

–Projection Utility (extension)

--Spatial Analyst (extension)

--Geoprocessing (extension)

--Grid Pig Tools (extension) see: http://webgis.wr.usgs.gov/gridpig.htm

--grid_cliptopoly.ave (Avenue script)

These and/or updated versions are available at ESRI’s website in the support|_downloads section, or at cost from ESRI.

http://www.esri.com/

specifically http://support.esri.com/index.cfm?fa=downloads.gateway

A geologic map of a ground water modelling region of interest, the Mahantango Creek Watershed of east-central Pennsylvania shows a paired anticline and syncline involving nine stratigraphic units:

II.1 Obtaining a geologic map and DEM and getting both into the same projection.

II.1.1 Obtaining and re-projecting a geologic map

The geologic map for the Mahantango site above came from the Pennsylvania SWAP dataset, several statewide coverages that have been digitized and georectified into the same Albers projection. Other geologic maps and coverages from other projections have been re-projected into the Albers projection.

After digitizing a desired coverage and noting the projection of the source map, the following routine in ArcView can be followed to get the coverage into a desired projection (in this example, from a UTM source to the desired Albers projection). Re-projecting coverages requires the Projection Utility of ArcView. The Projection Utility is added to ArcView as an extension by going to File|_ Extensions... and selecting it from the list of available extensions (if it is not available it must be purchased and/or downloaded, and installed in your ArcView EXT32 directory). The Projection Utility is then launched under File|_ArcView Projection Utility. Step 1 is to identify the shape file to be re-projected. Step 2 is to identify the input coordinate system. If a projection file exists for the shapefile (shapefile.prj), then the options here will be greyed out and the values from the prj file will be used. Step 3 is to identify the target output projection parameters. Albers (and a few other projections) must be handled as a custom projection. Then you select the datum and the projection as shown here:

Step 4 is to identify the output filename.

[Note1: I was able to re-project polygon (geology) shapefiles, but not point themes. I got an “error writing shapefile” message when attempting to re-project a point theme. This problem was solved for me by asking a colleague to re-project the point theme using ArcGIS (Arc8).]

[Note 2: It is unclear to me what happens with the greyed out values in step 2. I tried experimenting by projecting from UTM to Albers with a projection file and without (i.e. inserting the UTM input projection parameters manually) and got the same result, even though the “echo” of the parameter box just prior to processing showed the greyed out value (which in this case was defaulted to World Mercator due to an “unnamed” datum in the proj file). Therefore I am assuming a basic rule: if a prj file exists, it overrules all parameters entered manually.]

II.1.2 Obtaining, re-projecting, and resampling a digital elevation model (DEM).

Obtaining: A digital elevation model (DEM) for the Mahantango site was obtained from the USGS.

To obtain a DEM from the USGS seamless coverage, go to :

http://seamless.usgs.gov

and cut-out the desired DEM using their online viewer. The following are notes on how to use their viewer/DEM extractor:

Go to “View and order data sets – United States Viewer”.

[ Note: The screen resolution (viewing size) on your computer must be at least 1024 x 768 in order to see all options in the program, including the download option. To change your screen resolution, on your computer go to:

My Computer

|_Control Panel

|_Display

|_Settings

Select 1024 x 768 and hit O.K.. After screen goes blank and resets, click “Yes” within 15 seconds to avoid having the screen revert back to its original size. ]

In the viewer, use zoom tools to find your region of interest, then use the download tool to crop (i.e. square off) a DEM. Accept the defaults of NED (National Elevation Datum) and NLCD (National Land Coverage Datum) Raster downloads. The DEMs are on 30 meter lat-lon grid downloaded in a zip file. Unzip the zip file.

Re-projecting and importing the grid: The lat-lon grid of the original USGS DEM was projected into the Albers projection of the Mahantango project.

ArcView 3.2 does not have grid re-projecting and re-sampling routines. There may be a support extension to ArcView 3.2 to handle this; check the “downloads” section of ESRI. I had a colleague re-project these grids from the source projection (latitude-longitude, i.e. unprojected) to the Albers projection using ArcGIS (Arc8). The re-projected Albers grid was provided in Arc’s ASCII format.

To import an *.asc grid into ArcView, the Spatial Analyst extension must be on (selected under File|_extensions: if it is not available it must be purchased and/or downloaded, and installed in your ArcView EXT32 directory). Spatial Analyst allows import and export of ASCII grids (*.asc) under File|_Import/Export Data Source. Spatial Analyst also enables you to add a grid as a theme to a view. For View|_Add Theme, it adds “grid data source” to the default choices (“feature data source” and “image data source”) in the Data Source Types drop-down list of the Add Theme form. When adding a grid as a theme to a view, it reads the sta.adf file in a directory equal to the grid name.

Re-sampling: The original 30 meter resolution of the grid was considered too large, unnecessary for the scale of the regional Mahantango watershed analysis. The grid was resampled to a 90 meter resolution.

This underscores that copious data is available for defining DEMs at most project scales. To resample a grid to a smaller (90 meter) resolution, add the Grid Pig Tools to ArcView, including britecon.avx and gridpig.avx extensions, and g2i.dll to ArcView’s Bin32 directory. The extensions need to be activated in ArcView (selected under File|_extensions), which will add the “Grid Pig Tools” menu option to the menu bar. After selecting the grid to be resampled, select Grid Pig tools|_Grid Resample from the menu bar. The re-sampling routine will present an input box for the output cell size ( the current cell size of the grid is highlighted), and a pick arrow for the desired resampling method (with “nearest neighbor”, “bilinear interpolation”, and “cubic convolution” as choices). After clicking OK, the routine processes the new grid and imports it as a theme to the existing project.

II.2 Clipping the geologic map and DEM to a desired boundary (in this case a watershed boundary).

The geologic map above was clipped to the region of the Mahantango and Pine Creek combined watersheds. These watersheds were defined by processing the Digital Elevation Model (DEM) from the USGS with a hydrologic routing routine called TauDEM written by David Tarbotan (Tarboton, 1997).

TauDEM works with a DEM and a coincident (i.e. similarly projected) river polyline shapefile to define watershed boundaries. The program and directions for its use are available at:

http://hydrology.neng.usu.edu/taudem/

II.2.1 Clipping the geologic map theme to the watershed polygon theme

The geologic map theme for the Mahantango site was clipped to the boundary of the Mahantango watershed.

The Geoprocessing extension is required to clip (crop) a theme to the extents of another theme. Once the extension is activated (selected under File|_extensions: if it is not available it must be purchased and/or downloaded, and installed in your ArcView EXT32 directory), select the Geoprocessing wizard in the View menu. “Clip-one-theme-based-on-another” clips a theme to a border defined by another theme. When selected, a form is presented to specify 1) the input theme to clip, 2) the polygon overlay theme (the “clip to” theme), and 3) the output filename (a shape file). Use the watershed boundary shape file for step #2 to “clip” the geologic map theme to the extents of the watershed.

II.2.2 Clipping the DEM theme to a bounding box enclosing the watershed polygon theme

A bounding rectangle (polygon) theme was created around the Mahantango watershed, allowing extra room beyond the watershed in case changes to the boundary shape was required. The DEM grid was then clipped to the bounding rectangle.

To clip a grid to a polygon theme, you need the ArcView script (i.e. Avenue script) Grid_ClipToPoly.ave loaded into the project, compiled and run. Copy grid_cliptopoly.ave (the Avenue script) to ArcView’s scripts directory (typically: C:\ESRI\AV_GIS30\ARCVIEW\Samples\scripts). To compile and run the script in ArcView: 1) close the “view” and go to “scripts”, 2) select NEW button and generate script1, 3) Open script1 to get the “scripts” menus, 4) select Script|_Load System Script and load Grid.CliptoPoly, 5)select Script|_compile, 6) open the “view” and tile it with the “script” (not the project) window, 7) in the “view” window, select the grid theme you want to clip, 8) select Script|_Run from the menus. [Note: I initially got an “active document” error (i.e. it couldn’t find the themes to process because the “active document” wasn’t a view). You must make the “view” active, then click the “script” (not the project) window and run the script. If you are tiling: DO NOT TILE THE PROJECT WINDOW, ONLY THE VIEW AND SCRIPT WINDOWS, otherwise the script won’t work.] The script prompts, “Which polygon theme is the clipping theme?” Select the bounding theme from the pick-list. The script then prompt, “Clip grid on inside or outside of graphic?” I suggest using the outside of the graphic (i.e. the boundary polygon) to ensure the entire feature is within the grid. The script will then produce a new grid “clipped” to the bounding region.

The clipped DEM grid theme is then exported using File|_Export Data Source, selecting “ASCII raster” for the export file type, selecting the grid (from the list of available grids) to be exported, and specifying an output filename.

[Note: Though the script will clip a grid enclosing the entire graphic, it is usually a good idea to make a bounding rectangle around the graphic of interest and use the bounding rectangle to clip the grid.]

II.3 Converting the Geologic Map shapefile theme to a grid theme ( the integer geologic map raster known as geology-by-numbers).

The geologic map polygon theme for Mahantango was converted to a grid theme based on the stratigraphic “name” theme, however the resultant integers were not in stratigraphic order. This map was used for the projection example below. A similar grid map was produced, but using an ordered stratigraphic field for the conversion, resulting in stratigraphically ordered integers in the maps. This map was used for the folding example below.

Making a grid theme from a polygon theme requires that a theme (grid or polygon) of the desired extent exist (this can be the source theme itself). Though the cell and grid size can be set during the conversion, it is desirable to have a grid of the desired dimensions exist beforehand, particularly if you want to match the two grids, as in the case of a geology grid corresponding to a DEM grid. The geology shapefile theme is converted to a grid theme by selecting the geology theme and clicking on Theme|_Convert to Grid. The steps of this routine are 1) provide the output filename, 2) set both the “output grid extent” and “output grid cell size” to “same as [resampled DEM gridname]” by using the pick arrow and selecting the grid name of the resampled DEM. [Note: the cell size, number of columns, and number of rows should then automatically match the DEM grid], 3) select the conversion field to use for the integers (see next paragraph), and 4) answer “no” or “yes” (the answer is irrelevant) to the prompt, “join feature attributes to the grid?” The routine will generate a new grid theme, and then ask whether it should be added to the view (answer “yes” so that it can be selected and exported with File|_Export Data Source).

When ArcView converts a shapefile polygon theme to a raster grid theme, the integer values assigned to the grid theme for each polygon in the shapefile theme are based on the conversion export field, which is selected from a list of available fields in the database of the shapefile. Selecting “name” (the stratigraphic unit name) for this field works well (a unique integer is assigned for each stratigraphic unit), however the stratigraphic units are not in numerical order and the integers assigned to each stratigraphic name are essentially not user-controlled and thus not meaningful. If a new field is created in the shapefile database that supplies integers in stratigraphic order, and this field is used for step 3 above, it will be much easier to interpret the integers in the Viewxsxn program. To add a “strat” field, select the geologic map theme and go to Theme|_Table. With the table menus, select Table|_start editing. Then select Edit|_Add Field. Name the field (e.g. “Strat”), its type (number), and the number of decimal places ( 0 for integer). Click on the “selection” arrow-button. Using the shift key, select all the polygon records of a given stratigraphic name. Click on the “edit” I-beam/arrow button. Enter the stratigraphically-ordered layer number to represent this formation (set of polygons) by clicking the field’s (“strat’s) cell address for each selected record and placing the integer in each. For Viewxsxn, number the strata with “1" being the uppermost (usually youngest) layer and the total number of layers entered for the bottommost (usually oldest) layer. Deselect using the “deselect all” button. Repeat this process for each stratigraphic unit, then select Table|_stop editing. Save the edits by responding “yes” to the prompt. Close the table with File|_close, or with the upper right “x” on the table.

The geology grid theme is then exported using File|_Export Data Source, selecting “ASCII raster” for the export file type, selecting the grid (from the list of available grids) to be exported, and specifying an output filename.

II.4 Using Viewxsxn to generate a solid model

The geologic map and DEM for Mahantango were imported into Viewxsxn and folding and projection cross sections were generated.

Viewxsxn is started and the geology and DEM grids are read-in as instructed in section 2.0 introduction and 2.1. For a projection the Fold and Projection Parameters form was completed as follows:

In this case, the ArcView geology grid was exported converting on the “name” field, but because the stratigraphy was not ordered, the stratigraphic unit numbers have no relation to stratigraphic order. The stratigraphic unit numbers were related to the corresponding formation names only by visually scanning the output grid for each formation and jotting down the number. Similar integer maps might be constructed where a color number was exported for each formation and the color number is known, but needs to be stratigraphically ordered. In this case, the stratigraphy for the site ordered from top to bottom is listed with the order number (left) set equal to ( = ) the stratigraphic unit number “randomly” assigned by ArcView (right): Llewelyn (1=8), Pottsville (2=9), Mauch Chunk (3=6), Pocono (4=7), Spechty Kopf (5=5), Duncannon (6=1), Sherman Creek (7=3), Irish Valley (8=2), and Trimmers Rock (9=4). On the form, the order in which the stratigraphic unit numbers were encountered by Viewxsxn (which reads the map like a book) was 1,7,5,9,6,3,2,4,8. As a result, the assignments of the layer numbers were given to each of the stratigraphic unit numbers in the order they were encountered.

The user-completed form for each stratigraphic unit (left) has these top assignments (right):

1 6

7 4

5 5

9 2

6 3

3 7

2 8

4 9

8 1

After the stratigraphy is specified, clicking the arrow button causes Viewxsxn to calculate the average stratigraphic thicknesses based on outcrop widths. When these have been adjusted and/or accepted, the user selects either a folding or projection reconstruction with the radio button. In this case, the projection routine is selected:

When the projection is calculated, the dimensions of the resulting grids appear in the main menu:

The column row coordinates for a north-south cross section through column 139 from rows 1 to 209 were entered on the form, with equivalent coordinates (118936, 195552) to (118936, 176832) calculated by Viewxsxn based on the grid domain.

The cross section corresponds to the geologic map / cross section line:

When View was selected, the corresponding graph generated by Viewxsxn in Grapher shows the profile:

Notice that the dip is projected on both limbs of the anticline; the projection only depends upon the dip, the younging direction, and the distance to the surface contact. Many folds in the Valley and Ridge are chevron in nature with curvature localized to the hinge of the fold, and straight limbs. Projection works well for these types of folds. However, the depth of the projection is unrealistic: either the beds return to the surface forming a syncline, or they terminate into thrust faults at depth.

By trial and error, a set of parameters were determined to model an anticline-syncline fold pair entered onto the Fold and Projection Parameters form:

The x-offset (120200) and y-offset (185000) were determined by noting a coordinate on the limb of the fold midway between the anticlinal and synclinal axial traces, with the anticline in the positive y-direction and the syncline in the negative y-direction. This coordinate should be thought of as the origin for the mathematical generation of the fold. Though this fold shows the x-axis essentially east (N90°E), it is actually parallel to the strike (N85°E). The positive (toward the anticline) and negative (toward the syncline) y-axes should be thought of as perpendicular to strike. The z-offset (-1400), folding amplitude A (1700), and fold wavelength (20000) were determined by trial and error. The resulting fold cross section along the same line as the projection was generated with the View button:

A good technique to determine the three parameters is to initially set the z-offset too low, keeping it fixed while adjusting the amplitude and wavelength until they match the overlying outcrop connections. When the z-offset is too low, it is easy to find the contact on the surface, connected to the folding surface with vertical lines. An example with the same parameters as above, but with the z-offset set too low (-2400) is shown below:

Folds generated by stacking similar horizons [calculated by the same mathematical formula, but with different z-offsets] onto a basal horizon are similar folds: the stratigraphic thickness of any layer varies throughout the fold, thicker near the hinges and thinner on the limbs. In this case, the “depth thickness” is constant throughout while the stratigraphic thickness changes. Some folds in Pennsylvania are parallel folds, particularly on the Allegheny Plateau where anticlines and synclines are gentle warpings rather than tight, thrust-related folds. In future versions of Viewxsxn, I hope to incorporate faults and different (and better) fold shape models, including:

-parallel folds

-folds with wavelength/amplitude interdependencies

-spline fits to both bedding and horizon elevation data

-Fourier series generated fold shapes (see Stabler, 1968).

Note: If a geologic map is exported from ArcView based on the field of stratigraphically ordered integers, and read into Viewxsxn, the stratigraphic unit (left) and top numbers (right) are equal. Notice that the layer numbers are encountered in the file in the exact same order as with the “random” map. Figure shows form from preceding versions of Viewxsxn.

I currently have in legacy form a grid interpolator for generating grids from digitized, hand-drawn structure contour maps. I plan to incorporate this as a tool in Viewxsxn. I am also currently working on fold-model parameter optimization routines.