readme.txt

PREAMBLE: This folder contains the numerical simulations (Model Instances) and the software called FlowPaths (Model Program) reported in Inverse Method to Determine the Hydraulic Conductivity Parameter from a Velocity Field using Graph Theory by Michael E. Mont-Eton, Steffen Borgwardt, and David C. Mays. All software contained in this folder runs using Matlab R2021a for Windows. It is copyrighted by Michael Mont-Eton and available under the Creative Commons license Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0). The example simulation with two-dimensional matrix size 4x4 is as presented in the paper. Other simulations with matrix sizes of 4x4 through 16x16 are presented in a MATLAB workspace as described below. The tutorial provides instructions for performing various analyses of the data, and new data may be generated as described below.

CONTENTS of HydroShare FlowPaths_Resource folder

1. Files in the root directory:

readme.txt: (this file)
tutorial.pdf: Step-by-step tutorials for FlowPaths.
variables.pdf: Alphabetical listing of each variables name, units, description, class, and syntax.

2. Folders in the root directory:

/A.Documentation: License files for third-party MATLAB code used in FlowPaths.

/B.FlowPaths_Code_Files: The software programs (Model Program Aggregation) written using the MATLAB platform needed to run FlowPaths.

/C.MATLAB_Workspace_DataFiles: The MATLAB data files containing inputs and outputs of the FlowPaths simulations. This folder contains three subfolders as follows:
		
/C.1_4x4_Paper_Example_data: Input and output files for the 4x4 example.

/C.2_RecursiveStability_4x4_data: Files for recursive stability analysis of the 4x4 example.

/C.3_MultiFlow_Simulation_data: Input and output files for the multiple run simulation. This folder contains one subfolder as follows:

3. Files in each folder:

/A.Documentation

licols_license.txt : License file for the third-party MATLAB function licols.m.

rsgene2D_license.txt : License file for the third-party MATLAB function rsgene2D.m.

/B.FlowPaths_Code_Files

Note, see variables.pdf for definitions of the variables used below.

adjmat1b4.m
A subroutine used in DIRGRAPH4G4_Wght_MDST_LogErr_Z3.m that transforms the 2-dimensional m x n grid (rows x columns) of the initial flow matrix into an N x N graph adjacency matrix A, where N = mn is the number of cells in the grid, which is also equal to the number of vertices in the graph G. Each (i,j) entry in A is either 0 or 1, depending on whether the row entry (j) cell is within a Manhattan distance of one to the row (i) entry cell. Equivalently, if the jth entry is in the Von Neumann neighborhood (the 4 cardinal cells) of the ith entry then the (i,j) value is 1.

AllMaxZetaStat_v2.m 
A stand-alone program that compiles post-simulation analysis of the results from FlowPaths for each matrix size for the worst result of the Zeta error from all of the ten realizations for each Vdp number (i.e., each degree of heterogeneity) from the full data file of FlowPaths results. It accesses the tables of input and output Zeta data embedded in the multiple matrix size simulation set. The results are stored in a new MATLAB workspace and are also plotted as a function of Vdp versus Zeta for each Vdp number  in each matrix size. The purpose is to give the user an overall summary of the performance of FlowPaths and to identify simulations which had significantly lower accuracy. If the input variable verb2 is set to 1 (true), then a figure will be produced for each matrix size, showing the range of Zeta and ZetaMax across all of the realizations of increasing heterogeneity.

DIRGRAPH4G4_Wght_MDST_LogErr_Z3.m 
This program is the inverse solution generator. Its name derives from the purpose of creating the weighted digraph from the cells of the flow matrix grid, using a Minimum Directed Spanning Tree (MDST), with a log-error analysis of the result (Kfinal) compared to the initial K array.

DIRGRAPH4G4_Wght_MDST_LogErr_Z4.m
This program is slightly different from the Z3 version, and is intended to be used in the wrapper program WRPR_ISOLSEQ_MDST_v5_NP6_APZ_Test.m. It provides an output of the adjacency matrix, which can easily be stored in a single cell, whereas the former program automatically saves the Digraph and the adjacency matrix in a mat-file. 

InVarTabl_Mkr_4RecStats_SinglSol_Chk_v2.m 
This program organizes the variables needed to run the recursive stability checker (WRPR_FlowPaths_SnglVdp_Stats_v4.m) into a MATLAB table named InVarTable, stored in a MATLAB workspace file.

licols.m 
This subroutine performs a reduction of the over-determined array of q* coefficients (the combination of specific discharges entering and leaving each vertex along any path) produced by DIRGRAPH4G4_Wght_MDST_LogErr_Z3.m to a linearly independent matrix of q* coefficients with full rank. Its name is based on Linearly Independent Columns. The license for this third-party code is provided in licols_license.txt as mentioned above.

LoadDefault_FlowPaths_Wrapper_Invars.m
This program loads the variables needed to create a new set of hydraulic conductivities to be used in the FlowPaths system. All plot and save variables are turned off (verbose=0; XpP=0; Xportall=0).

LoadDefault_FlowPaths_Wrapper_Invars_4Plot.m
This program is the same as above but changes the plot and save variables to verbose=1; XpP=1; Xportall=1.

LoadDefault_Multi_Matrix_Vdp_Invars.m 
This subroutine loads the default variables needed to populate the workspace for creating the multi-matrix-size MATLAB data structure (see below). It is called by MakeBlankMasterStrucDataFile.m. It creates a workspace file that contains the default variables, which is opened by MakeBlankMasterStrucDataFile.m.

LoadDefault_Single_Sim_AnyUID_Invars.m 
This program loads the default settings for running a check of any FlowPaths single simulation, based on the UID that the user chooses to analyze. The output includes the initial hydraulic conductivity array along with boundary conditions, and all the variables needed to run the forward and inverse solver contained in the wrapper program WRPR_SINGLESOLCHECK_MDST_APZ.m.

MakeBlankMasterStrucDataFile_v2.m 
This program creates the initial MATLAB data structure M_InvSolSeq_DataFile, which is needed to run a multi-matrix-size, multi-Vdp, multi-realization FlowPaths simulation. Its name derives from its purpose to make a blank master structure data file. The default variables loaded into the workspace from LoadDefault_Multi_Matrix_FlowPath_Invars.m and the user-designated string variables TLF, CF, clx, and cly are needed as inputs.

MakeSingle_InvSolSeq_DataFile.m
This program creates an initial data structure for a single flow matrix: M_InvSolSeq_Single_DataFile. It is used to store the data from a single FlowPaths simulation. It is necessary to run the VDP_HK_gen_2D_MM7_R6.m program first to generate the inputs for this program, which are the independent variables Ncol, Nrow, and CF. The structure HKGEN from the output of VDP_HK_gen_2D_MM7_R6.m is also needed as input. Please see variables.pdf for definitions. The output of this file is slightly different from the MakeBlankMasterStrucDataFile.m program in that there is only one row of fields and there is only one UID.

PBN_MDST_sub_v2.m 
This program is the path-finding subroutine used in the FlowPaths inverse solver DIRGRAPH4G4_Wght_MDST_LogErr_Z3.m. The input variables of note are: DGMAT1, src, and snk. DGMAT1 is the digraph that is developed from the specific discharges through the whole flow matrix. src is the source vertex chosen from the high-head boundary vertices. snk is the target, or sink vertex chosen from the low-head boundary vertices.  The program defines the MDST between the source vertex (s) and the sink vertex (t), as the spanning tree, rooted at s, choosing the minimum-weight edge (as a function of the specific discharges going through them) at each possible branching. The program then defines two types of paths through the MDST sub-digraph: the chord paths and the one non-chord path. The first set is a group of paths that are each concatenated from three parts: (1) a path through the MDST from the source vertex src to the vertex (u) at the end of a branch, (2) the chord with the smallest weight from (u) to its nearest neighbor vertex (v), and (3) the least-weighted path from (v) to snk. The latter path is the only path that traverses the MDST from src to snk without going through any chords. The output of the program, which is passed back to DIRGRAPH4G4_Wght_MDST_LogErr_Z3.m, is a set of vertices for each path from s to t.

PR_FWD_qGEN_4p_v3_Z5.m 
This program solves the forward groundwater equation problem, using the initial K array (i.e., hydraulic conductivity array) and the boundary head conditions. It produces the block-centered head and the face-centered specific discharge fields. The specific discharge field is broken out in terms of the left faces, the right faces, the top faces, and the bottom faces, which is similar to how MODFLOW organizes the specific discharges in its forward solution. The solution method is based on an alternating direction implicit algorithm, adapted from the work of Esfandiari (2017, p. 429, https://doi.org/10.1201/9781315152417). Its name is derived from Peaceman-Rachford Alternating Direction Implicit Method. The specific discharge arrays (qL, qR, qT, and qB) are the data used by the inverse solver, FlowPaths.

rsgene2D.m 
This program generates two-dimensional spatially-correlated random fields of hydraulic conductivity with a normal distribution and an exponential auto-covariance. These random fields are interpreted as the log of the hydraulic conductivity, rendering lognormally distributed random fields of hydraulic conductivity. Its name derives from Random Surface generator in 2D. The correlation lengths (isotropic) are set at 0.05 cm. The license for this third-party code is provided in rsgene2D_license.txt as mentioned above.

VDP_HK_gen_2D_MM7_R7.m 
This program generates the initial hydraulic conductivity field in units of m/day, using rsgene2D.m with the Vdp variable (the Reservoir Heterogeneity Index) determining the degree of heterogeneity. The Gaussian distribution of heterogeneity is centered on the default absolute permeability of 0.00000001 square centimeters. The programs name derives from using Vdp (a number from zero to one, where zero represents absolute uniformity of the hydraulic conductivity field and one represents infinite heterogeneity) to generate the HK (hydraulic conductivity) variable (an m by m array of heterogeneous isotropic hydraulic conductivities, where m is the number of rows in the array) in 2D, along with the version descriptors as a suffix. The suggested range of Vdp is from 0.01 to 0.98. For details, please see tutorial.pdf.

VDP_HK_gen_2D_Wrapper_R6.m 
This program generates a large set of HK arrays based on the lower and upper criteria for the Vdp, having ten random realizations for each heterogeneity value of Vdp. It calls the VDP_HK_gen_2D_MM7_R7.m program as a subroutine. Its name derives from its usage as a wrapper function- one that contains a series of subroutine calls to generate iterations of the main subroutine. 

WRPR_FlowPaths_SnglVdp_Stats_v5.m
This statistical analysis program determines whether the inverse solution is stable over several recursions of executing the forward and inverse solvers (PR_FWD_qGEN_4p_v3_Z5.m and DIRGRAPH4G4_Wght_MDST_LogErr_Z3.m, respectively). The relative accuracy between the inverse solution and the initial K data, defined as Zeta, is used as the independent variable in the statistical analysis routine from the MATLAB library called cusumtest.m, which produces 5 outputs including the series of recursive residuals, which is a set of independent variables (e.g., Galpin and Hawkins, 1984). Its name derives from the use of a single Vdp in the FlowPaths system that produces a statistical inference at the 95% confidence level for whether the particular solution is stable, and thus meets one of the requirements for the problem being well-posed.

WRPR_FlowPaths_SnglVdp_Stats_v6.m
This program is a variant of the one above, in that the Zeta values are log-shifted back from the ones used in version 5. We provide it for doing an additional analysis of the cumulative sums of squares test for stability. 

WRPR_ISOLSEQ_MDST_v5_NP6_APZ.m 
This program is a wrapper function that computes results several matrix sizes, one at a time, the full range of Vdp values and the K-fields, the forward solver, and the inverse solver to produce solutions for over 10,000 flow matrices. All data is stored in a mat-file named MSolSetFileName.mat that contains a structure of structures named M_InvSolSeq_DataFile. Its name derives from Inverse Solution Sequence using the MDST. It is provided for information only, as users will want to use the _Test version described next.

WRPR_ISOLSEQ_MDST_v5_NP6_APZ_Test.m
This program is the same as WRPR_ISOLSEQ_MDST_v5_NP6_APZ.m, except that the names of the mat-file and the structure have been changed (MySolSetFileName.mat and My_M_InvSolSeq_DataFile, respectively), so that a reviewer may check the results of the original simulations without overwriting them. For details, please see tutorial.pdf.

WRPR_SINGLESOL_Plot_MDST_APZ.m
This program controls both the forward and inverse solvers (PR_FWD_qGEN_4p_v3_Z5.m, and DIRGRAPH4G4_Wght_MDST_LogErr_Z3.m), with full plot and save functionality.

WRPR_SINGLESOLCHECK_MDST_APZ.m 
This program is also a wrapper function, but it only runs through one sequence of forward and inverse solutions, using a single K field array as input. The output is placed into a data structure inside a workspace file. Its name derives from just doing a single solution check, using the MDST approach in the inverse solver subroutine. For details, please see tutorial.pdf.

ZetaTest_SingleMat.m
This program compares the users test of all of the entries (cells) of hydraulic conductivity to the data in the original mat-file for a full set of realizations for one flow matrix size. If all the differences between Zeta statistics are zero, it prints Zeta Equal. If they are not all zero, it prints Zeta Not Equal.

/C.MATLAB_Workspace_DataFiles

/C.1_4x4_Paper_Example_data

WS_InvSol_InitVars_4x4_2023-03-12_Vdp_0.5_Rz_01_DToday_20230611.mat
This MATLAB workspace file (otherwise known as a mat-file) contains the initial variables (ds, HB, HK, HL, HR, HT, itermax, tol, verbose, XpP) needed to run the program WRPR_SINGLESOLCHECK_MDST_APZ.m, except for the current folder variable, CF, which tells the computer where to store the results and figures (if any) generated by its execution. The file name is WS  for workspace; InvSol for inverse solution; InitVars for initial variables; 4x4_2023-03-12_Vdp_0.5_Rz_01, which is the unique identifier for the 4x4 example given in the paper; and Dtoday_20230611  for the date that the initial dataset was created.

WS_InvSol_OutVars_4x4_2023-03-12_Vdp_0.5_Rz_01_DToday_20230611.mat
This MATLAB workspace file contains the results of the forward solver and the inverse solver for the 4x4 example matrix, as called by WRPR_SINGLESOLCHECK_MDST_APZ.m. It also contains all of the initial variables, including CF. The naming convention differs from the initial variable mat-file by the use of OutVars.

/C.2_RecursiveStability_4x4_data

DefRecStabInvars.mat
This mat-file contains the default stability input variables needed to run the stability checking program (WRPR_FlowPaths_SnglVdp_Stats_v4.m).

WKspace_DGmat_Digraph-4x6_20230710.mat
This mat-file is created during the execution of DIRGRAPH4G4_Wght_MDST_LogErr_Z3.m, which occurs inside the WRPR_FlowPaths_SnglVdp_Stats_v4.m program. It contains the digraph of the flow matrix (DGmat1) and the weighted adjacency matrix (A1), as the specific discharges going from vertex u to vertex v in the digraph. WKspace_DGmat_ specifies that the workspace contains information about the digraph for the 4x6 full flow matrix, and the last string is the date it was created. Since the code overwrites this mat-file every time it repeats, this mat-file is the last one produced by the recursive stability checker.

WS_CUSUM_StatAnalysis_OutVars_UID_4x4_2023-03-12_Vdp_0.5_Rz_01_Iter_50_DToday_20230710.mat
This mat-file contains all of the results of the MATLAB cumulative sum test for stability for the 4x4 example. Most notably, the result of the test was that NullHyp = 0, which means that the simulation passed the stability test. WS_CUSUM_StatAnalysis_OutVars indicates that this is the results workspace for the CUSUM test for stability. The rest of the strings are the same as for the input variables workspace.

WS_InVars_StatAnalysis_UID_(4x4_2023-03-12_Vdp_0.5_Rz_01)_Iter_50_DToday_20230710.mat
This MATLAB workspace file contains the initial variables and the InVarTable table (created by InVarTabl_Mkr_4RecStats_Chk.m) for the analysis done by the program WRPR_FlowPaths_SnglVdp_Stats_v4.m, of the 4x4 example used in the paper. The string WS_InVars_StatAnalysis says the mat-file contains the initial variables for the statistical analysis for stability; UID_(4x4_2023-03-12_Vdp_0.5_Rz_01) is the UID descriptor; Iter_50 is the number of recursions; and the last string is the date of the organization of the workspace.

/C.3_MultiFlow_Simulation_data

FlowPaths_UID_Tables_from_MasterDataFile_DToday_20230623.mat
This mat-file contains two structures: (1) The master data file, M_InvSolSeq_DataFile, for the 13 flow matrix sizes (4x4 to 16x16), over a range of Vdp values from 0.01 to 0.98, having ten random realizations each. (2) The structure, UIDTableStruct, which has two fields: CMS and UIDTable. Each row has one entry for the Core Matrix Size and one entry for the MATLAB table containing 980 entries (98 Vdp values times 10 realizations) and 13 headings. The first column of the table is the UID, which enables a faster search of the results than using M_InvSolSeq_DataFile. The string FlowPaths_UID_Tables_from_MasterDataFile indicates that this mat-file was generated from the master data file structure.

Master_InvSolSeq_Datafile_4to16_2023-03-19.mat
This mat-file also contains the master data file structure M_InvSolSeq_DataFile. The string Master_InvSolSeq_Datafile_ indicates that it is the overall data set for the multi-flow-matrix set of FlowPaths simulations; 4to16 indicates that the square flow matrices start with 4 rows and continue on to 16 rows; 2023-03-19 is the date that the multi-flow-matrix simulation was initiated.

My_Default_InvSolSeq_Invars.mat
A set of default inputs for kick-starting the FlowPaths wrapper for the forward and inverse solutions.

REFERENCES

Esfandiari, R. (2017). Numerical Methods for Engineers and Scientists Using MATLAB (2nd ed.). Boca Raton, Florida, USA: Taylor and Francis.

Galpin, J. S., & Hawkins, D. M. (1984). The Use of Recursive Residuals in Checking Model Fit in Linear-Regression. American Statistician, 38(2), 94-105. doi:10.2307/2683242