GMS:T-PROGS: Difference between revisions

From XMS Wiki
Jump to navigationJump to search
No edit summary
 
(32 intermediate revisions by 4 users not shown)
Line 1: Line 1:
{{Infobox Model |
{{Infobox Model |
|name=          T-PROGS
|name=          T-PROGS
|screenshot=    Image:Matset.gif
|screenshot=    Image:Matset.png
|model_type=    transition probability geostatistics on borehole data  
|model_type=    transition probability geostatistics on borehole data  
|developer=    Steven Carle  
|developer=    Steven Carle  
Line 9: Line 9:
{{T-PROGS Links}}
{{T-PROGS Links}}


GMS includes an interface to the T-PROGS software developed by Steven Carle.  The T-PROGS software is used to perform transition probability geostatistics on borehole data.  The output of the T-PROGS software is a set of N [[GMS:Material Set|material sets]] on a 3D grid.  Each of the material sets is conditioned to the borehole data and the materials proportions and transitions between the boreholes follows the trends observed in the borehole data.  These material sets can be used for [[GMS:Stochastic Modeling#Stochastic Applications|stochastic simulations]] with MODFLOW.  A sample material set generated by the TPROGS software is shown below.  The T-PROGS software can also be used to generate multiple input data sets for the [[GMS:HUF Package|HUF package]].
GMS includes an interface to the T-PROGS (Transition PRObability GeoStatistical) software developed by Steven Carle.  The T-PROGS software is used to perform transition probability geostatistics on borehole data.  The output of the T-PROGS software is a set of N [[GMS:Material Set|material sets]] on a 3D grid.  Each of the material sets is conditioned to the borehole data and the materials proportions and transitions between the boreholes following the trends observed in the borehole data.  These material sets can be used for [[GMS:Stochastic Modeling#Stochastic Applications|stochastic simulations]] with MODFLOW.  A sample material set generated by the T-PROGS software is shown below.  This software can also be used to generate multiple input datasets for the [[GMS:HUF Package|HUF package]].
 
The T-PROGS model can be added to a [http://www.aquaveo.com/software/gms-pricing paid edition] of GMS.
[[Category:Link to Store]]


==T-PROGS Interface==
==T-PROGS Interface==


{|
{|
| [[Image:pathlines.png|frame|none|''MODFLOW solution for a T-PROGS grid. Black lines are head contours and blue lines are pathlines for an extraction well.''|250px]]
|-valign="top"
| [[Image:markov_graphs.png|frame|none|''Observed transition probabilities (circles) and Markov chain model (solid lines) for a set of borehole data with four materials.''|250px]]
| [[Image:pathlines.png|thumb|none|350px|MODFLOW solution for a T-PROGS grid. Black lines are head contours and blue lines are pathlines for an extraction well.]]
| [[Image:markov_graphs.png|thumb|none|350px|Observed transition probabilities (circles) and Markov chain model (solid lines) for a set of borehole data with four materials.]]
|}
|}


The T-PROGS software utilizes a transition probability-based geostatistical approach to model spatial variability by 3-D Markov Chains, set up indicator co-kriging equations, and formulate the objective function for simulated annealing.  The transition probability approach has several advantages over traditional indicator kriging methods. First, the transition probability approach considers asymmetric juxtapositional tendencies, such as fining-upwards sequences.  Second, the transition probability approach has a conceptual framework for incorporating geologic interpretations into the development of cross-correlated spatial variability.  Furthermore, the transition probability approach does not exclusively rely on empirical curve fitting to develop the indicator (cross-) variogram model.  This is advantageous because geologic data are typically only adequate to develop a model of spatial variability in the vertical direction.  The transition probability approach provides a conceptual framework to geologic insight into a simple and compact mathematical model, the Markov chain.  This is accomplished by linking fundamental observable attributes mean lengths, material proportions, anisotropy, and juxtapositioning with Markov chain model parameters.
The T-PROGS software utilizes a transition probability-based geostatistical approach to model spatial variability by 3-D Markov Chains, set up indicator co-kriging equations, and formulate the objective function for simulated annealing.  The transition probability approach has several advantages over traditional indicator kriging methods. First, the transition probability approach considers asymmetric juxtapositional tendencies, such as fining-upwards sequences.  Second, the transition probability approach has a conceptual framework for incorporating geologic interpretations into the development of cross-correlated spatial variability.  Furthermore, the transition probability approach does not exclusively rely on empirical curve fitting to develop the indicator (cross-) variogram model.  This is advantageous because geologic data are typically only adequate to develop a model of spatial variability in the vertical direction.  The transition probability approach provides a conceptual framework to geologic insight into a simple and compact mathematical model, the Markov chain.  This is accomplished by linking fundamental observable attributes—mean lengths, material proportions, anisotropy, and juxtapositioning – with Markov chain model parameters.


The first step in using T-PROGS is to import a set of borehole data.  The borehole data are then passed to a utility within T-PROGS called GAMEAS that computes a set of [[GMS:Vertical Markov Chain|transition probability curves]] as a function of lag distance for each category for a given sampling interval.  A sample set of measured transition probability curves are shown by the dashed lines in the following figure.
The first step in using T-PROGS is to import a set of borehole data.  The borehole data are then passed to a utility within T-PROGS called GAMEAS that computes a set of [[GMS:Vertical Markov Chain|transition probability curves]] as a function of lag distance for each category for a given sampling interval.  A sample set of measured transition probability curves are shown by the dashed lines in the following figure.


Each curve represents the transition probability from material j to material k.  The transition probability t<sub>jk</sub>(h) is defined by:
Each curve represents the transition probability from material j to material k.  The transition probability ''t<sub>jk</sub>(h)'' is defined by:


<!--[[Image:t-progs_eq1.png]]
<!--[[Image:t-progs_eq1.png]]
<math>t_{jk}(\mathbf{h})= \text{Pr}\{k \text{ occurs at } \mathbf{x} + \mathbf{h}\ |\ j \text{ occurs at } \mathbf{x}\}</math>-->
<math>t_{jk}(\mathbf{h})= \text{Pr}\{k \text{ occurs at } \mathbf{x} + \mathbf{h}\ |\ j \text{ occurs at } \mathbf{x}\}</math>-->
[[image:tprogseq1.jpg]]<div style="float: right;"><ref>
:[[image:tprogseq1.jpg]]<div style="float: right;"><ref>
{{citation  
{{citation  
|last=Carle
|last=Carle
Line 37: Line 41:
}}</ref></div>
}}</ref></div>


where x is a spatial location, h is the lag (separation vector), and j,k denote materials.  Note that the curves on the diagonal represent auto-transition probabilities, and the curves on the off-diagonal represent cross-transition probabilities.
where ''x'' is a spatial location, ''h'' is the lag (separation vector), and ''j,k'' denote materials.  Note that the curves on the diagonal represent auto-transition probabilities, and the curves on the off-diagonal represent cross-transition probabilities.


The next step in the analysis is to develop a Markov Chain model for the vertical direction that fits the observed vertical transition probability data.  The Markov Chain curves are shown as solid lines in the preceeding figure.  Mathematically, a Markov chain model applied to one-dimensional categorical data in a direction Φ assumes a matrix exponential form:
The next step in the analysis is to develop a Markov Chain model for the vertical direction that fits the observed vertical transition probability data.  The Markov Chain curves are shown as solid lines in the preceeding figure.  Mathematically, a Markov chain model applied to one-dimensional categorical data in a direction Φ assumes a matrix exponential form:
Line 43: Line 47:
<!-- [[Image:t-progs_eq2.png]]
<!-- [[Image:t-progs_eq2.png]]
<math>\mathbf{T}(h_\phi)=\text{exp}(\mathbf{R}_{\phi}h_\phi)</math>-->
<math>\mathbf{T}(h_\phi)=\text{exp}(\mathbf{R}_{\phi}h_\phi)</math>-->
[[image:tprogseq2.jpg]]<div style="float: right;"><ref>
:[[image:tprogseq2.jpg]]<div style="float: right;"><ref>
{{citation  
{{citation  
|last=Carle
|last=Carle
Line 54: Line 58:
}}</ref></div>
}}</ref></div>


where h denotes a lag in the direction Φ, and RΦ denotes a transition rate matrix
where ''h'' denotes a lag in the direction ''Φ'', and '''' denotes a transition rate matrix


[[Image:rate_matrix.gif]]
:[[Image:rate_matrix.png]]


with entries r<sub>jk</sub>,f representing the rate of change from category j to category k (conditional to the presence of j) per unit length in the direction Φ.  The transition rates are adjusted to ensure a good fit between the Markov Chain model and the observed transition probability data.
with entries ''r<sub>jk</sub>,f'' representing the rate of change from category ''j'' to category ''k'' (conditional to the presence of ''j'') per unit length in the direction ''Φ''.  The transition rates are adjusted to ensure a good fit between the Markov Chain model and the observed transition probability data.


Once the Markov chain is developed for the z direction from the borehole data, a model of spatial variability must be developed for the x and y directions.  Borehole data are typically not sufficiently dense in these directions.  However, the [[GMS:Strike Dip Markov Chain|x and y-direction Markov chains]] can be developed by assuming that the juxtapositional tendencies and the proportions observed in the vertical direction also hold true in the horizontal directions.  The modeler then provides an estimate of the ratio of the mean lengths in the x and y directions relative to the z direction, and the transition rate matrices for the x and y directions can be formulated.  The x, y, and z Markov chains are converted into a continuous 3D Markov chain using the MCMOD utility within T-PROGS.
Once the Markov chain is developed for the z direction from the borehole data, a model of spatial variability must be developed for the x and y directions.  Borehole data are typically not sufficiently dense in these directions.  However, the [[GMS:Strike Dip Markov Chain|x and y-direction Markov chains]] can be developed by assuming that the juxtapositional tendencies and the proportions observed in the vertical direction also hold true in the horizontal directions.  The modeler then provides an estimate of the ratio of the mean lengths in the x and y directions relative to the z direction, and the transition rate matrices for the x and y directions can be formulated.  The x, y, and z Markov chains are converted into a continuous 3D Markov chain using the MCMOD utility within T-PROGS.
Line 64: Line 68:
In the final phase of setting up a transition probability analysis using T-PROGS, the modeler creates a grid, specifies the number of model instances (N), and launches the TSIM utility.  The TSIM code uses the 3D Markov chain to formulate both indicator cokriging equations and an objective function for simulated annealing.  It generates [[GMS:Stochastic Modeling#Stochastic Applications|stochastic simulations]] using a combination of modified versions of the GSLIB codes SISIM and ANNEAL.
In the final phase of setting up a transition probability analysis using T-PROGS, the modeler creates a grid, specifies the number of model instances (N), and launches the TSIM utility.  The TSIM code uses the 3D Markov chain to formulate both indicator cokriging equations and an objective function for simulated annealing.  It generates [[GMS:Stochastic Modeling#Stochastic Applications|stochastic simulations]] using a combination of modified versions of the GSLIB codes SISIM and ANNEAL.


==T-PROGS OPTIONS==
TPROGS allows roughly up to 3.5 million cells to be included in a simulation. Model instabilities may appear if more that 3.5 million cells are used.
''T-PROGS Options dialog with and without borehole data.''
 
==T-PROGS Materials==
{|
{|
| [[Image:t-progs_wo_boreholes.png|frame|none|''Without boreholes'']]
| [[Image:t-progs_wo_boreholes.png|thumb|none|300 px|''T-PROGS Materials'' dialog without borehole data.]]
| [[Image:t-progs_w_boreholes.png|frame|none|''With boreholes'']]
| [[Image:t-progs_w_boreholes.png|thumb|none|300 px|''T-PROGS Materials'' dialog with borehole data.]]
|}
|}


When a user selects the '''''New Simulation''''' command to initialize a T-PROGS simulation, the ''T-PROGS Options'' dialog appears. If boreholes do not exist in the model, an unconditioned simulation will be generated. In this case, the user selects the materials to be used and a corresponding background material. The '''[[GMS:Materials|Materials Editor]]''' button enables users to quickly create, delete, or rename materials and the material list is automatically regenerated.  The upper part of the dialog lists the materials in the boreholes.  The first column of toggles indicates which materials are to be used in the analysis.  By default, all materials associated with the boreholes are selected.  These toggles are necessary since it is possible that there may be materials defined in the materials list that are not associated with boreholes.  The second column in the top section of the dialog lists the background material.  By default, the material type that had the predominant occurrence in the boreholes (greatest proportion) is marked as the background material.  When defining the transition probability data in the next section, the input parameters do not need to be edited for the background material.  The parameters for this material are automatically adjusted to balance the equations.
When having selected the '''New Simulation''' command to initialize a T-PROGS simulation, the ''T-PROGS Boreholes'' dialog appears. Here select to use all boreholes or only the boreholes in a particular folder.
 
===HGUs===
Boreholes use materials to define both soils and [[GMS:Borehole_Hydrogeologic_Units|HGUs]]. HGUs can be used to group several materials into one hydrogeologic unit. The ''T-PROGS Materials'' dialog lets choosing to use the materials in all HGUs or those from just one HGU. If ''Use all HGUs'' is selected, T-PROGS will use the HGU information on the boreholes and ignore the soil information. If ''Use one HGU'' is selected, T-PROGS will use the soil information on the boreholes for the soils that are in the selected HGU. This feature can be used to limit the portion of the boreholes that are used in the T-PROGS simulation.
 
===Material List===
If boreholes do not exist in the model, an unconditioned simulation will be generated. In this case, in the ''T-PROGS Materials'' dialog select the materials to be used and a corresponding background material. The upper part of the dialog lists the materials in the boreholes.  The first column of toggles indicates which materials are to be used in the analysis.  By default, all materials associated with the boreholes are selected.  These toggles are necessary since it is possible that there may be materials defined in the materials list that are not associated with boreholes.  The second column in the top section of the dialog lists the background material.  By default, the material type that had the predominant occurrence in the boreholes (greatest proportion) is marked as the background material.  When defining the transition probability data in the next section, the input parameters do not need to be edited for the background material.  The parameters for this material are automatically adjusted to balance the equations.


===Background Material===
===Background Material===
Line 78: Line 89:


===Azimuth===
===Azimuth===
''Azimuths for T-PROGS and 3D Grid.''
{|multiple image
{|multiple image
| [[Image:azimuth1.png|frame|none|''T-Progs''|150px]]
| [[Image:azimuth1.png|frame|none|Azimuths for T-PROGS.|150px]]
| [[Image:azimuth2.png|frame|none|''3D Grid''|175px]]
| [[Image:azimuth2.png|frame|none|Azimuths for 3D Grid.|175px]]
|}
|}


The user also enters an azimuth in this dialog. The azimuth determines the orientation of the primary directions of the depositional trends in the strike/dip directions.  These trends generally are aligned with the primary directions of horizontal flow in the aquifer.  Theoretically, the azimuth can be oriented independently from the grid orientation.  However, in practice, if the grid and azimuth orientations are offset by more than about 40<sup>o</sup>, checkerboard patterns appear in the indicator array results.  Hence, the azimuth orientation is set equal to the grid orientation by default.  However, the grid angle is defined counterclockwise, and the azimuth angle is clockwise.  Therefore, if the grid angle is 40<sup>o</sup>, then the azimuth angle will be –40<sup>o</sup> by default.  If there is anisotropy in the xy plane, the azimuth angle should be set to the principle direction of the anisotropy.  If anisotropy is not present, this angle should be coincident with the x-axis (the rows or j-direction) of the grid.
Also enter an azimuth in this dialog. The azimuth determines the orientation of the primary directions of the depositional trends in the strike/dip directions.  These trends generally are aligned with the primary directions of horizontal flow in the aquifer.  Theoretically, the azimuth can be oriented independently from the grid orientation.  However, in practice, if the grid and azimuth orientations are offset by more than about 40<sup>o</sup>, checkerboard patterns appear in the indicator array results.  Hence, the azimuth orientation is set equal to the grid orientation by default.  However, the grid angle is defined counterclockwise, and the azimuth angle is clockwise.  Therefore, if the grid angle is 40<sup>o</sup>, then the azimuth angle will be –40<sup>o</sup> by default.  If there is anisotropy in the xy plane, the azimuth angle should be set to the principle direction of the anisotropy.  If anisotropy is not present, this angle should be coincident with the x-axis (the rows or j-direction) of the grid.


===Material Limit===
===Material Limit===
Line 91: Line 101:


==Generating Material Sets with T-PROGS==
==Generating Material Sets with T-PROGS==
The underlying equations solved by the T-PROGS software require an orthogonal grid with constant cell dimensions (X, Y, and Z).  The delta X values can be different from the delta Y and delta Z values, and the delta Y values can be different from the delta Z values, but all cells must have the same change in X, Y, and Z dimensions. The MODFLOW model is capable of using the [[GMS:LPF Package|Layer Property Flow (LPF) Package]] with the Material ID option for assigning aquifer properties.  With this option, each cell in the grid is assigned a material id and the aquifer properties (Kh, Kv, etc.) associated with each material are automatically assigned to the layer data arrays for the LPF package when the MODFLOW files are saved. The T-PROGS software generates multiple material sets (arrays of material ids), each of which represents a different realization of the aquifer heterogeneity.  When running a MODFLOW simulation in [[GMS:Stochastic Modeling#Stochastic Applications|stochastic mode]], GMS automatically loads each of the N material sets generated by the T-PROGS software and saves N different sets of MODFLOW input files.  The N solutions resulting from these simulations can be read into GMS and used to perform [[GMS:Risk_Analysis_Wizard|risk analyses]] such as probabilistic capture zone delineation.
The underlying equations solved by the T-PROGS software require an orthogonal grid with constant cell dimensions (X, Y, and Z).  The delta X values can be different from the delta Y and delta Z values, and the delta Y values can be different from the delta Z values, but all cells must have the same change in X, Y, and Z dimensions. The MODFLOW model is capable of using the [[GMS:LPF Package|Layer Property Flow (LPF) package]] with the ''Material ID'' option for assigning aquifer properties.  With this option, each cell in the grid is assigned a material ID and the aquifer properties (Kh, Kv, etc.) associated with each material are automatically assigned to the layer data arrays for the LPF package when the MODFLOW files are saved. The T-PROGS software generates multiple material sets (arrays of material IDs), each of which represents a different realization of the aquifer heterogeneity.  When running a MODFLOW simulation in [[GMS:Stochastic Modeling#Stochastic Applications|stochastic mode]], GMS automatically loads each of the N material sets generated by the T-PROGS software and saves N different sets of MODFLOW input files.  The N solutions resulting from these simulations can be read into GMS and used to perform [[GMS:Risk_Analysis_Wizard|risk analyses]] such as probabilistic capture zone delineation.


=== One-Layer MODFLOW Grids ===
=== One-Layer MODFLOW Grids ===
[[Image:matset_samp.png|right|thumb|''Sample material set generated by the T-PROGS software.  Blue lines represent particle path lines and black lines are head contours.''|275px]]
[[Image:matset_samp.png|right|thumb|Sample material set generated by the T-PROGS software.  Blue lines represent particle path lines and black lines are head contours.|275px]]


Although MODFLOW is a three-dimensional model, a majority of the MODFLOW models constructed by typical users are 2D models consisting of one model layer.  There are several reasons why 2D models are so common.  One reason is that many of these models are regional models where the aquifer thickness is very small compared to the lateral extent of the model.  As a result, the flow directions are primarily horizontal and little improvement is gained by adding multiple layers to the model.  Even with local scale models, the aquifer thickness is often small enough that one-layer models are considered adequate.  2D models are also attractive due to the simplicity of the model increased computational efficiency.  One of the problems associated with using multiple layers for MODFLOW models with unconfined aquifers is that as the water table fluctuates, the upper cells may go dry.  These cells will not rewet even if the water table subsequently rises, unless the rewetting option has been selected in the flow package (BCF, LPF, or HUF).  The rewetting issues can often be avoided with a one-layer model.
Although MODFLOW is a three-dimensional model, a majority of the MODFLOW models constructed are 2D models consisting of one model layer.  There are several reasons why 2D models are so common.  One reason is that many of these models are regional models where the aquifer thickness is very small compared to the lateral extent of the model.  As a result, the flow directions are primarily horizontal and little improvement is gained by adding multiple layers to the model.  Even with local scale models, the aquifer thickness is often small enough that one-layer models are considered adequate.  2D models are also attractive due to the simplicity of the model increased computational efficiency.  One of the problems associated with using multiple layers for MODFLOW models with unconfined aquifers is that as the water table fluctuates, the upper cells may go dry.  These cells will not rewet even if the water table subsequently rises, unless the rewetting option has been selected in the flow package (BCF, LPF, or HUF).  The rewetting issues can often be avoided with a one-layer model.


When developing a one-layer model, the modeler must determine how to distribute the hydraulic conductivity values within the layer.  One option is to assume a homogenous aquifer; this is typically a gross over-simplification since aquifers are usually highly heterogeneous.  Therefore, a common approach is to delineate zones of hydraulic conductivity by examining the subsurface stratigraphic data.  In many cases, these data are in the form of borehole logs.  These borehole logs often exhibit substantial heterogeneity and don’t always exhibit definitive trends between adjacent boreholes.  Furthermore, the boreholes are often clustered with large regions of the model lacking any borehole data.  The modeler then faces a difficult task of trying to determine a rational approach to delineating two-dimensional zones of hydraulic conductivity based on complex 3D borehole data.
When developing a one-layer model, the modeler must determine how to distribute the hydraulic conductivity values within the layer.  One option is to assume a homogenous aquifer; this is typically a gross over-simplification since aquifers are usually highly heterogeneous.  Therefore, a common approach is to delineate zones of hydraulic conductivity by examining the subsurface stratigraphic data.  In many cases, these data are in the form of borehole logs.  These borehole logs often exhibit substantial heterogeneity and don’t always exhibit definitive trends between adjacent boreholes.  Furthermore, the boreholes are often clustered with large regions of the model lacking any borehole data.  The modeler then faces a difficult task of trying to determine a rational approach to delineating two-dimensional zones of hydraulic conductivity based on complex 3D borehole data.
Line 106: Line 116:
Using transition probability geostatistics with MODFLOW models results in two basic limitations.  First, the underlying stochastic algorithms used by the T-PROGS software are formulated such that the MODFLOW grid must have uniform row, column, and layer widths.  The row width can be different from the column width, but each row must have the same width.  This results in a uniform orthogonal grid.  While MODFLOW grids are orthogonal in x and y, the layer thickness is allowed to vary on a cell-by-cell basis.  This makes it possible for the layer boundaries to accurately model the ground surface and the tops and bottoms of aquifer units.  If a purely orthogonal grid is used, irregular internal and external layer boundaries must be simulated in a stair-step fashion either by varying material properties or by activating/inactivating cells via the IBOUND array.  A second limitation is that in order to get a high level of detail in the simulated heterogeneity, the grid cell dimensions are generally kept quite small.  This can result in difficulties in the vertical dimension.  The large number of layers with small layer thicknesses near the top of the model generally ensures that many of the cells in this region will be at or above the computed water table elevation (for simulations involving unconfined aquifers).  As a result, these cells will undergo many of the numerical instabilities and increased computational effort issues associated with cell wetting and drying.
Using transition probability geostatistics with MODFLOW models results in two basic limitations.  First, the underlying stochastic algorithms used by the T-PROGS software are formulated such that the MODFLOW grid must have uniform row, column, and layer widths.  The row width can be different from the column width, but each row must have the same width.  This results in a uniform orthogonal grid.  While MODFLOW grids are orthogonal in x and y, the layer thickness is allowed to vary on a cell-by-cell basis.  This makes it possible for the layer boundaries to accurately model the ground surface and the tops and bottoms of aquifer units.  If a purely orthogonal grid is used, irregular internal and external layer boundaries must be simulated in a stair-step fashion either by varying material properties or by activating/inactivating cells via the IBOUND array.  A second limitation is that in order to get a high level of detail in the simulated heterogeneity, the grid cell dimensions are generally kept quite small.  This can result in difficulties in the vertical dimension.  The large number of layers with small layer thicknesses near the top of the model generally ensures that many of the cells in this region will be at or above the computed water table elevation (for simulations involving unconfined aquifers).  As a result, these cells will undergo many of the numerical instabilities and increased computational effort issues associated with cell wetting and drying.


[[Image:huf_sample.gif|right|frame|''Sample HUF stratigraphy data generated by the T-PROGS software.  Blue lines represent particle path lines and water table.''|300px]]
[[Image:huf_sample.png|thumb|right|475 px|Sample HUF stratigraphy data generated by the T-PROGS software.  Blue lines represent particle path lines and water table.]]


The Hydrogeologic Unit Flow (HUF) package  released with MODFLOW 2000 makes it possible to overcome both of these limitations resulting in a powerful mechanism for incorporating transition probability geostatistics in MODFLOW simulations.  With the HUF package, the modeler is allowed to input the vertical component of the stratigraphy in a grid-independent fashion.  The stratigraphy data are defined using a set of elevation and thickness arrays.  The first array defines the top elevation of the model.  The remaining arrays define the thicknesses of a series of hydrogeologic units, starting at the top and progressing to the bottom of the model.  For each array of thicknesses, many of the entries in the array may be zero.  This makes it possible to simulate complex heterogeneity, including pinchouts and embedded lenses that would be difficult to simulate with the LPF and BCF packages.
The Hydrogeologic Unit Flow (HUF) package  released with MODFLOW 2000 makes it possible to overcome both of these limitations resulting in a powerful mechanism for incorporating transition probability geostatistics in MODFLOW simulations.  With the HUF package, the modeler is allowed to input the vertical component of the stratigraphy in a grid-independent fashion.  The stratigraphy data are defined using a set of elevation and thickness arrays.  The first array defines the top elevation of the model.  The remaining arrays define the thicknesses of a series of hydrogeologic units, starting at the top and progressing to the bottom of the model.  For each array of thicknesses, many of the entries in the array may be zero.  This makes it possible to simulate complex heterogeneity, including pinchouts and embedded lenses that would be difficult to simulate with the LPF and BCF packages.


The T-PROGS interface in GMS includes an option for integrating transition probability geostatistics results with the HUF package.  The basic approach used by the option is to overlay a dense background grid on the MODFLOW grid and run T-PROGS on the background grid.  A set of HUF arrays is then extracted from the background grid for use with the MODFLOW model.  To use this option, user should first create a MODFLOW grid with the desired number of layers and the layer elevations should be interpolated to match the aquifer boundaries.  The row and column widths are uniform but the layer thicknesses may vary from cell to cell.  Then, when TSIM is launched, the HUF option should be selected.  GMS then generates a background grid that encompasses the MODFLOW grid.  The rows and columns of this grid match the MODFLOW grid but the layer thicknesses are uniform and relatively thin, resulting in a much greater number of layers than the MODFLOW grid.  The user specifies the number of layers in this background grid.  A T-PROGS simulation is then performed to get a set of material sets on the background grid.  Each of the material sets in the T-PROGS output is then transferred from the background grid to a set of HUF elevation/thickness arrays.  The HUF top elevation array is set equal to the top of the MODFLOW grid.  The thickness arrays are then found by searching through the background grid to find the bottom elevations of contiguous groups of indicators.  The elevations from these groups are then added to an appropriate elevation array in the HUF input.  The resulting set of HUF input arrays are listed in GMS [[GMS:The GMS Screen|Project Explorer]].  By clicking on each item in the [[GMS:Project Explorer|Project Explorer]], the selected set of HUF arrays are loaded into the HUF package and the corresponding stratigraphy is displayed in the GMS window.  The multiple HUF input arrays can be used to perform a [[GMS:Stochastic Modeling#Stochastic Applications|stochastic simulation]].
The T-PROGS interface in GMS includes an option for integrating transition probability geostatistics results with the HUF package.  The basic approach used by the option is to overlay a dense background grid on the MODFLOW grid and run T-PROGS on the background grid.  A set of HUF arrays is then extracted from the background grid for use with the MODFLOW model.  To use this option, first create a MODFLOW grid with the desired number of layers and the layer elevations should be interpolated to match the aquifer boundaries.  The row and column widths are uniform but the layer thicknesses may vary from cell to cell.  Then, when TSIM is launched, the HUF option should be selected.  GMS then generates a background grid that encompasses the MODFLOW grid.  The rows and columns of this grid match the MODFLOW grid but the layer thicknesses are uniform and relatively thin, resulting in a much greater number of layers than the MODFLOW grid.  Specify the number of layers in this background grid.  A T-PROGS simulation is then performed to get a set of material sets on the background grid.  Each of the material sets in the T-PROGS output is then transferred from the background grid to a set of HUF elevation/thickness arrays.  The HUF top elevation array is set equal to the top of the MODFLOW grid.  The thickness arrays are then found by searching through the background grid to find the bottom elevations of contiguous groups of indicators.  The elevations from these groups are then added to an appropriate elevation array in the HUF input.  The resulting set of HUF input arrays are listed in GMS [[GMS:The GMS Window|Project Explorer]].  By clicking on each item in the [[GMS:Project Explorer|Project Explorer]], the selected set of HUF arrays are loaded into the HUF package and the corresponding stratigraphy is displayed in the GMS window.  The multiple HUF input arrays can be used to perform a [[GMS:Stochastic Modeling#Stochastic Applications|stochastic simulation]].


==Notes==
==Notes==
Line 119: Line 129:
[[Category:T-PROGS]]
[[Category:T-PROGS]]
[[Category:Stochastic]]
[[Category:Stochastic]]
[[Category:T-PROGS Dialogs]]
[[Category:GMS Stochastic Modeling Dialogs]]
[[Category:Equations|T-Progs]]
[[Category:External Links]]

Latest revision as of 17:17, 29 May 2020

T-PROGS
T-PROGS
T-PROGS Screenshot
Model Info
Model type transition probability geostatistics on borehole data
Developer Steven Carle
Documentation T-Progs Manual
Tutorials T-PROGS Tutorials

T-PROGS
Links
T-PROGS
Strike Dip Markov Chain
Vertical Markov Chain
TSIM

GMS includes an interface to the T-PROGS (Transition PRObability GeoStatistical) software developed by Steven Carle. The T-PROGS software is used to perform transition probability geostatistics on borehole data. The output of the T-PROGS software is a set of N material sets on a 3D grid. Each of the material sets is conditioned to the borehole data and the materials proportions and transitions between the boreholes following the trends observed in the borehole data. These material sets can be used for stochastic simulations with MODFLOW. A sample material set generated by the T-PROGS software is shown below. This software can also be used to generate multiple input datasets for the HUF package.

The T-PROGS model can be added to a paid edition of GMS.

T-PROGS Interface

MODFLOW solution for a T-PROGS grid. Black lines are head contours and blue lines are pathlines for an extraction well.
Observed transition probabilities (circles) and Markov chain model (solid lines) for a set of borehole data with four materials.

The T-PROGS software utilizes a transition probability-based geostatistical approach to model spatial variability by 3-D Markov Chains, set up indicator co-kriging equations, and formulate the objective function for simulated annealing. The transition probability approach has several advantages over traditional indicator kriging methods. First, the transition probability approach considers asymmetric juxtapositional tendencies, such as fining-upwards sequences. Second, the transition probability approach has a conceptual framework for incorporating geologic interpretations into the development of cross-correlated spatial variability. Furthermore, the transition probability approach does not exclusively rely on empirical curve fitting to develop the indicator (cross-) variogram model. This is advantageous because geologic data are typically only adequate to develop a model of spatial variability in the vertical direction. The transition probability approach provides a conceptual framework to geologic insight into a simple and compact mathematical model, the Markov chain. This is accomplished by linking fundamental observable attributes—mean lengths, material proportions, anisotropy, and juxtapositioning – with Markov chain model parameters.

The first step in using T-PROGS is to import a set of borehole data. The borehole data are then passed to a utility within T-PROGS called GAMEAS that computes a set of transition probability curves as a function of lag distance for each category for a given sampling interval. A sample set of measured transition probability curves are shown by the dashed lines in the following figure.

Each curve represents the transition probability from material j to material k. The transition probability tjk(h) is defined by:

Tprogseq1.jpg

where x is a spatial location, h is the lag (separation vector), and j,k denote materials. Note that the curves on the diagonal represent auto-transition probabilities, and the curves on the off-diagonal represent cross-transition probabilities.

The next step in the analysis is to develop a Markov Chain model for the vertical direction that fits the observed vertical transition probability data. The Markov Chain curves are shown as solid lines in the preceeding figure. Mathematically, a Markov chain model applied to one-dimensional categorical data in a direction Φ assumes a matrix exponential form:

Tprogseq2.jpg

where h denotes a lag in the direction Φ, and denotes a transition rate matrix

Rate matrix.png

with entries rjk,f representing the rate of change from category j to category k (conditional to the presence of j) per unit length in the direction Φ. The transition rates are adjusted to ensure a good fit between the Markov Chain model and the observed transition probability data.

Once the Markov chain is developed for the z direction from the borehole data, a model of spatial variability must be developed for the x and y directions. Borehole data are typically not sufficiently dense in these directions. However, the x and y-direction Markov chains can be developed by assuming that the juxtapositional tendencies and the proportions observed in the vertical direction also hold true in the horizontal directions. The modeler then provides an estimate of the ratio of the mean lengths in the x and y directions relative to the z direction, and the transition rate matrices for the x and y directions can be formulated. The x, y, and z Markov chains are converted into a continuous 3D Markov chain using the MCMOD utility within T-PROGS.

In the final phase of setting up a transition probability analysis using T-PROGS, the modeler creates a grid, specifies the number of model instances (N), and launches the TSIM utility. The TSIM code uses the 3D Markov chain to formulate both indicator cokriging equations and an objective function for simulated annealing. It generates stochastic simulations using a combination of modified versions of the GSLIB codes SISIM and ANNEAL.

TPROGS allows roughly up to 3.5 million cells to be included in a simulation. Model instabilities may appear if more that 3.5 million cells are used.

T-PROGS Materials

T-PROGS Materials dialog without borehole data.
T-PROGS Materials dialog with borehole data.

When having selected the New Simulation command to initialize a T-PROGS simulation, the T-PROGS Boreholes dialog appears. Here select to use all boreholes or only the boreholes in a particular folder.

HGUs

Boreholes use materials to define both soils and HGUs. HGUs can be used to group several materials into one hydrogeologic unit. The T-PROGS Materials dialog lets choosing to use the materials in all HGUs or those from just one HGU. If Use all HGUs is selected, T-PROGS will use the HGU information on the boreholes and ignore the soil information. If Use one HGU is selected, T-PROGS will use the soil information on the boreholes for the soils that are in the selected HGU. This feature can be used to limit the portion of the boreholes that are used in the T-PROGS simulation.

Material List

If boreholes do not exist in the model, an unconditioned simulation will be generated. In this case, in the T-PROGS Materials dialog select the materials to be used and a corresponding background material. The upper part of the dialog lists the materials in the boreholes. The first column of toggles indicates which materials are to be used in the analysis. By default, all materials associated with the boreholes are selected. These toggles are necessary since it is possible that there may be materials defined in the materials list that are not associated with boreholes. The second column in the top section of the dialog lists the background material. By default, the material type that had the predominant occurrence in the boreholes (greatest proportion) is marked as the background material. When defining the transition probability data in the next section, the input parameters do not need to be edited for the background material. The parameters for this material are automatically adjusted to balance the equations.

Background Material

Application of the transition probability approach involves the designation of a background material. The probabilistic constraints of the Markov chains make it unnecessary to quantify data for one category. Not only is it unnecessary, but it is futile to do so because values will be overwritten in order to satisfy constraints. Conceptually, the background material can be described as the material that “fills” in the remaining areas not occupied by other units. For example, in a fluvial depositional system, a floodplain unit would tend to occupy area not filled with higher-energy depositional units and would therefore be a logical choice for the background material.

Azimuth

Azimuths for T-PROGS.
Azimuths for 3D Grid.

Also enter an azimuth in this dialog. The azimuth determines the orientation of the primary directions of the depositional trends in the strike/dip directions. These trends generally are aligned with the primary directions of horizontal flow in the aquifer. Theoretically, the azimuth can be oriented independently from the grid orientation. However, in practice, if the grid and azimuth orientations are offset by more than about 40o, checkerboard patterns appear in the indicator array results. Hence, the azimuth orientation is set equal to the grid orientation by default. However, the grid angle is defined counterclockwise, and the azimuth angle is clockwise. Therefore, if the grid angle is 40o, then the azimuth angle will be –40o by default. If there is anisotropy in the xy plane, the azimuth angle should be set to the principle direction of the anisotropy. If anisotropy is not present, this angle should be coincident with the x-axis (the rows or j-direction) of the grid.

Material Limit

One limitation for both the cases with and without boreholes is that a maximum of five materials can be used in the T-PROGS algorithm. This limitation was imposed to keep the data processing and user-interface reasonably simple. Although five materials present a limitation, borehole data can generally be easily condensed down to five or fewer materials. Furthermore since this is a stochastic approach, which is based on probability, the detail generated with numerous materials is rarely justifiable anyway. In addition, as the number of materials increase, the ratio of process time to detail becomes inefficient.

Generating Material Sets with T-PROGS

The underlying equations solved by the T-PROGS software require an orthogonal grid with constant cell dimensions (X, Y, and Z). The delta X values can be different from the delta Y and delta Z values, and the delta Y values can be different from the delta Z values, but all cells must have the same change in X, Y, and Z dimensions. The MODFLOW model is capable of using the Layer Property Flow (LPF) package with the Material ID option for assigning aquifer properties. With this option, each cell in the grid is assigned a material ID and the aquifer properties (Kh, Kv, etc.) associated with each material are automatically assigned to the layer data arrays for the LPF package when the MODFLOW files are saved. The T-PROGS software generates multiple material sets (arrays of material IDs), each of which represents a different realization of the aquifer heterogeneity. When running a MODFLOW simulation in stochastic mode, GMS automatically loads each of the N material sets generated by the T-PROGS software and saves N different sets of MODFLOW input files. The N solutions resulting from these simulations can be read into GMS and used to perform risk analyses such as probabilistic capture zone delineation.

One-Layer MODFLOW Grids

Sample material set generated by the T-PROGS software. Blue lines represent particle path lines and black lines are head contours.

Although MODFLOW is a three-dimensional model, a majority of the MODFLOW models constructed are 2D models consisting of one model layer. There are several reasons why 2D models are so common. One reason is that many of these models are regional models where the aquifer thickness is very small compared to the lateral extent of the model. As a result, the flow directions are primarily horizontal and little improvement is gained by adding multiple layers to the model. Even with local scale models, the aquifer thickness is often small enough that one-layer models are considered adequate. 2D models are also attractive due to the simplicity of the model increased computational efficiency. One of the problems associated with using multiple layers for MODFLOW models with unconfined aquifers is that as the water table fluctuates, the upper cells may go dry. These cells will not rewet even if the water table subsequently rises, unless the rewetting option has been selected in the flow package (BCF, LPF, or HUF). The rewetting issues can often be avoided with a one-layer model.

When developing a one-layer model, the modeler must determine how to distribute the hydraulic conductivity values within the layer. One option is to assume a homogenous aquifer; this is typically a gross over-simplification since aquifers are usually highly heterogeneous. Therefore, a common approach is to delineate zones of hydraulic conductivity by examining the subsurface stratigraphic data. In many cases, these data are in the form of borehole logs. These borehole logs often exhibit substantial heterogeneity and don’t always exhibit definitive trends between adjacent boreholes. Furthermore, the boreholes are often clustered with large regions of the model lacking any borehole data. The modeler then faces a difficult task of trying to determine a rational approach to delineating two-dimensional zones of hydraulic conductivity based on complex 3D borehole data.

As part of this research, we developed a technique for developing 2D zones of hydraulic conductivity from borehole logs using transition probability geostatistics. The technique is simple, fast, and preserves proportions and trends exhibited by the borehole data. The algorithm parses through each borehole and computes a predominant material at each borehole. When T-PROGS runs, the predominant material for each borehole is assigned to its corresponding location in the one-layer grid, and during the quenching process, simulations are conditioned to those data points.

Generating HUF Data with T-PROGS

Using transition probability geostatistics with MODFLOW models results in two basic limitations. First, the underlying stochastic algorithms used by the T-PROGS software are formulated such that the MODFLOW grid must have uniform row, column, and layer widths. The row width can be different from the column width, but each row must have the same width. This results in a uniform orthogonal grid. While MODFLOW grids are orthogonal in x and y, the layer thickness is allowed to vary on a cell-by-cell basis. This makes it possible for the layer boundaries to accurately model the ground surface and the tops and bottoms of aquifer units. If a purely orthogonal grid is used, irregular internal and external layer boundaries must be simulated in a stair-step fashion either by varying material properties or by activating/inactivating cells via the IBOUND array. A second limitation is that in order to get a high level of detail in the simulated heterogeneity, the grid cell dimensions are generally kept quite small. This can result in difficulties in the vertical dimension. The large number of layers with small layer thicknesses near the top of the model generally ensures that many of the cells in this region will be at or above the computed water table elevation (for simulations involving unconfined aquifers). As a result, these cells will undergo many of the numerical instabilities and increased computational effort issues associated with cell wetting and drying.

Sample HUF stratigraphy data generated by the T-PROGS software. Blue lines represent particle path lines and water table.

The Hydrogeologic Unit Flow (HUF) package released with MODFLOW 2000 makes it possible to overcome both of these limitations resulting in a powerful mechanism for incorporating transition probability geostatistics in MODFLOW simulations. With the HUF package, the modeler is allowed to input the vertical component of the stratigraphy in a grid-independent fashion. The stratigraphy data are defined using a set of elevation and thickness arrays. The first array defines the top elevation of the model. The remaining arrays define the thicknesses of a series of hydrogeologic units, starting at the top and progressing to the bottom of the model. For each array of thicknesses, many of the entries in the array may be zero. This makes it possible to simulate complex heterogeneity, including pinchouts and embedded lenses that would be difficult to simulate with the LPF and BCF packages.

The T-PROGS interface in GMS includes an option for integrating transition probability geostatistics results with the HUF package. The basic approach used by the option is to overlay a dense background grid on the MODFLOW grid and run T-PROGS on the background grid. A set of HUF arrays is then extracted from the background grid for use with the MODFLOW model. To use this option, first create a MODFLOW grid with the desired number of layers and the layer elevations should be interpolated to match the aquifer boundaries. The row and column widths are uniform but the layer thicknesses may vary from cell to cell. Then, when TSIM is launched, the HUF option should be selected. GMS then generates a background grid that encompasses the MODFLOW grid. The rows and columns of this grid match the MODFLOW grid but the layer thicknesses are uniform and relatively thin, resulting in a much greater number of layers than the MODFLOW grid. Specify the number of layers in this background grid. A T-PROGS simulation is then performed to get a set of material sets on the background grid. Each of the material sets in the T-PROGS output is then transferred from the background grid to a set of HUF elevation/thickness arrays. The HUF top elevation array is set equal to the top of the MODFLOW grid. The thickness arrays are then found by searching through the background grid to find the bottom elevations of contiguous groups of indicators. The elevations from these groups are then added to an appropriate elevation array in the HUF input. The resulting set of HUF input arrays are listed in GMS Project Explorer. By clicking on each item in the Project Explorer, the selected set of HUF arrays are loaded into the HUF package and the corresponding stratigraphy is displayed in the GMS window. The multiple HUF input arrays can be used to perform a stochastic simulation.

Notes