Home

User`s Guide to Parssim1 - The Institute for Computational

1. lt real L gt lt real L gt zMin zMax Minimal and maximal z points The minimal and maximal z points of the computational grid 31 5 4 2 Nonuniform Rectangular Grid Use the following when gridType is rectangular For each of the nx 1 x grid points give the following real L gt xp The x grid points j 1 For each of the ny 1 y grid points give the following real L gt yp The y grid points For each of the nz 1 z grid points give the following real L gt zp The z grid points 5 4 3 An XY rectangular Z variable Grid Use the following when gridType is xy rectangular For each of the nx 1 x grid points give the following real L gt xp The x grid points For each of the ny 1 y grid points give the following real L gt yp The y grid points j For each of the nx 1 ny 1 nz 1 z grid points give the following lt real of a data block L gt zp The z grid points See 84 7 3 for information on a data block 5 4 4 A Fully Logically Rectangular Grid Use the following when gridType is This indicates the beginning of a data block For each of the nx 1 ny 1 nz 1 x grid points give real of a data block L gt xp The x grid points See 84 7 3 for information on a data block 32 j 1 For each of the nx 1 ny 1 nz 1 y grid points give the following real of a data block L gt yp The y grid points See 84 7
2. zWellIndexLo Lower z index of well If wellType selects a reinjection well specify the lower z index number of the extraction part lt omit gt lt integer gt zWellIndexHi Higher z index of well If wellType selects a reinjection well specify the upper z index number of the extraction part 45 5 14 1 1 Flow SSS Flow This sectional unit is needed only if computeFlow 0 Bottomfole Producer Use if wellType selects a bottomhole producer real M L T 2 wellPres Well pressure potential The pressure potential of the well see 2 1 on page 6 for the meaning of potential Other Active Well Types The following lines appear only if wellType does not select a bottomhole producer lt integer gt nWellInterp f Number of interpolation times The number of interpolation times used in the definition of the wells as a piecewise continuous linear function of time Constant extrapolation is used before and after the first and last times thus a single interpolation time gives a constant condition The times must not decrease however two pairs with the same time represent a discontinuous jump 1 For each of the nWellInterp f interpolation times give the following real T gt welllnj f time Well time The time at which the next interpolation value is to be used real L 3 T gt wellInj f value Well injection or extraction rate The next interpolation value volume of injected or extracted fluid per u
3. For each of the nBCInterp_t interpolation times give the following real T bc t time BC time real 1 L 3 bc t value BC conc value S WELLS lt integer gt nWells Number of wells lt integer gt maxnWellInterp_f Maximum number of flow interpolation times lt integer gt maxnWellInterp t Maximum number of transport interpolation times SS Well Include one such sectional unit per well as specified by nWells 0111112113114 wellType Well type If the well is not inactive include the following information lt real L gt wellRadius Well radius lt integer gt lt integer gt hWellIndex The x and y indices of well lt integer gt zWellIndexLo Lower z index of well lt integer gt zWelllndexHi Higher z index of well omit integer zWellIndexLo Lower z index of well omit integer zWellIndexHi Higher z index of well SSS Flow This sectional unit is needed only if computeFlow 0 Use if wellType selects a bottomhole producer real M L T 2 1 gt wellPres Well pressure potential The following lines appear only if wellType does not select a bottomhole producer lt integer gt nWellInterp_f Number of interpolation times For each of the nWelllnterp f interpolation times give the following real T wellInj_f_time Well time real L 3 T gt wellInj_f_value Well injection or extraction rate SSS Transport This sectional unit
4. Gibbs free energy 5 11 13 15 Godunov Method 5 9 10 26 gravitational constant 6 gravity 6 33 GRID 31 grid 5 22 27 28 31 32 45 56 grid array 22 22 33 35 41 43 53 grid 000 56 gridType 31 31 32 H Component 37 H20 Component 37 half saturation constants 39 halfLife 40 halfSatConst 14 39 head 43 49 49 height 31 Henry s Law 6 8 37 hessFlag 29 Higher Order Godunov Method see Godunov Method histIndexHi t x 57 histIndexHi t x y z 90 histIndexHi t y 57 histIndexHi t z 57 histIndexLo t x 57 histIndexLo t x y z 90 histIndexLo t y 57 histIndexLo t z 57 History 50 histSpecie t 50 histType t 50 HOG see Godunov Method hWellIndex 45 hydraulic conductivity 23 hydraulicHead 43 49 49 ideal 36 ideal 36 immobile 8 33 37 inactive well 45 45 incompressible 5 7 36 86 index range 22 23 45 47 50 57 infile 18 25 inflow condition 43 44 INITIAL CONDITIONS 41 initial0ut 48 injection well 8 45 interface problem 7 27 49 interior point algorithm 5 11 17 29 ion exchange 11 iWorkspace_c 30 k 6 Kelvin 21 key word 19 21 52 53 BOUNDARY CONDITIONS 42 branching 40 40 by 34 34 cells 34 Chain 40 CHEMICAL SPECIES 37 Chemistry 28 36 57 CMM 28 Component 37 components 34 Conductivity 33 diagonal 33 34 Dispersion 35 face 33 35 Flow 27 30 41 43 46 49 GENERAL INFO 25 GRID 31 H Co
5. i e the La norm of the residuals the KKT conditions 18 must be less than or equal to the user specified tolerance This condition has different implications depending on the choice of equilibrium module however and we discuss here the SF and the RNSF Using the SF mass balance is exactly satisfied and the residuals measure deviation from op timality and complementarity as explained in 25 It is easy to see that in a system in which all species have concentrations which are sufficiently large satisfaction of the stopping criterion will approximately guarantee that for the ith reaction AG chmAbsTol 2 18 where G is the Gibbs free energy These results can be expressed as CNot i 57 exp chmAbsTol z 1 chmAbsTol 2 19 Ki j ey or equivalently IN s leno i Ki ITA Cei Ki TG ej s lt chmAbsTol 2 20 15 which clearly shows the role of chmAbsTol as a kind of relative error in this application For species present at small concentrations this bound will deteriorate In the RNSF on the other hand mass balance is only approximately satisfied by a given iterate and the absolute value of the mass balance error for any component is bounded by chmAbsTol The implication is that even a small specified tolerance could result in a large relative error for a component which is present in trace amounts only Note that in the RNSF the residuals consist entirely of component mass balances in the
6. A value of zero indicates that no sorption takes place A negative value indicates an immobile specie and the value of ajo is assumed to be 1 integer omit sorpType Sorption type If needed phaseDist gt 0 the sorption capacity number from 1 to nSorpTypes specified earlier This uses the order of the sorption capacities induced above SSS Chemistry This sectional unit is needed only if computeChemistry gt 0 Components and products are described differently The command line argument c can be used to see how a given data set is interpreted 37 5 8 1 1 Component Chemistry Specification Use this syntax if the subsection heading is Component H Component or H20 Component lt integer gt phaseldentity Participating phase The phase in which the component participates based on the organization of phases given in elsewhere lt real 1 gt compCharge Charge The intrinsic charge of the component Used to compute the intrinsic species charges and it can be valuable for debugging the stoichiometry 5 8 1 2 Product Chemistry Specification Use this syntax if the subsection heading is Product lt integer gt phaseldentity Participating phase The phase in which the product participates based on the organization of phases given in elsewhere For each of the nComps components give the following lt real 1 gt stoich Stoichiometry formula number The formula number for this product species
7. END COMMENT BEGIN COMMENT Causes input to be ignored until the matching END COMMENT command is read such commands are nested or the end of file is reached END COMMENT Terminates a BEGIN COMMENT command IGNORE LINES Follow by proper white space and an integer Causes that number of lines i e up to newline to be ignored from the input stream The line containing the integer is included in the count of lines ignored and information is ignored from that point on However the include command only is not ignored This command can be used to ignore preambles in included data files 19 IGNORE_WORDS Follow by proper white space and an integer Causes that number of words to be ignored from the input stream However the include command only is not ignored This command can be used to ignore preambles in included data files INCLUDE Follow by proper white space and a new file name maximum of 50 characters This causes subsequent reading to take place in the new file When the end of file is reached reading commences in the original file While there is no limit on the number of files that may be included there is a maximum level to which they may be nested LITERAL Terminates reading white space and other symbols This command is used to allow a special symbol at the beginning of a data string such as punctuation and comment symbol and command symbol This command is not ne
8. It is expected that several species may use the same o distribution and that often o If linear sorption is used and a species is mobile the concentration variable conc refers to the concentration in the fluid phase The total concentration T is then Ti aici ci so that o is essentially the Henry s Law constant An immobile species has T ajoicj Linear and nonlinear sorption can be simulated by the general chemistry routines If both the linear sorption model described here and the chemistry routines are used care must be taken in setting up the physical problem Unlike the general chemistry routines this linear sorption model does not treat the adsorbed or absorbed specie as a distinct specie It is recommended that only one of these two models be used at a time 2 3 3 Numerical Solution The advection and diffusion dispersion subproblems are solved independently using time splitting Three options can be selected for the solution of the advection subproblem 2 3 3 1 Basic Time Splitting Algorithm The algorithm used can be described by the following steps We let m denote the current time level and At the current time step size Some of the following sub problems may not be computed for some or all of the species as appropriate Advection We solve the equation Oc aioi as V GU S qei For the Characteristics Mixed Method we solve for T ajo c according to OT u VULT ai
9. could be well before this point in the file s Error on processor lt number gt An unspecified error occurred on the given processor Error opening input file lt file name gt The given file was not able to be opened Perhaps it is misspelled or there is an error in the directory tree path to the file Expression has nonnested parentheses A units expression cannot have nonnested parentheses Grid array expected A grid array see 4 7 3 on page 22 was expected but not found Perhaps a begin list symbol was left out Key not found lt list of words gt The key word from the given list was not found It is likely missing misspelled or out of order Premature beginning of SECTION command found A new section was begun before the last one was finished Premature beginning of SUBSECTION command found A new subsection was begun before the last one was finished 53 Premature beginning of SUBSUBSECTION command found A new sub subsection was begun before the last one was finished Premature end of data The data file ended before all appropriate data was read and processed Section key not found lt list of words gt The section name from the given list was not found It is likely missing misspelled or out of order Sorry there is an error in the computer code Refer the problem to someone who can debug the code This error message occurs generally only for errors that are relatively straightforward to resolve Sub
10. phaseIdentity Participating phase 1 compCharge Charge SS Component 2 B specieName Specie name 100 0 molecularWeight Molecular weight 0 0 phaseDist The alpha parameter in Henry s Law SSS Chemistry 1 phaseIdentity Participating phase 1 compCharge Charge SS Component 3 76 C specieName Specie name 200 0 molecularWeight Molecular weight 0 0 phaseDist The alpha parameter in Henry s Law SSS Chemistry 1 phaseIdentity Participating phase 1 compCharge Charge SS Component 4 D specieName Specie name 200 0 molecularWeight Molecular weight 0 0 phaseDist The alpha parameter in Henry s Law SSS Chemistry 1 phaseIdentity Participating phase 1 compCharge Charge SS Product 1 AB s specieName Specie name 200 molecularWeight Molecular weight 99 phaseDist The alpha parameter in Henry s Law SSS Chemistry 2 phaseIdentity Participating phase 1 0 1 0 0 0 0 0 st stoich 0 Equilibrium 0 0 pK Log 10 equilibrium constant SS Product 2 AC s specieName Specie name 300 molecularWeight Molecular weight 99 phaseDist The alpha parameter in Henry s Law SSS Chemistry 3 phaseIdentity Participating phase 1 0 0 0 1 0 0 0 stoich 0 Equilibrium 3 0103e 1 pK Log 10 equilibrium constant SS Product 3 DB s specieName Specie name 300 molecularWeight Molecular weight 99 phaseDist The alpha paramete
11. 15 26 27 29 38 39 43 45 48 55 u 18 18 21 54 v 18 18 ol 48 21 21 19 20 19 20 lt gt 25 22 22 Eau 190 C J 21 52 54 CC 54 20 20 BEGIN_COMMENT 19 19 20 END_COMMENT 19 19 20 54 IGNORE_LINES 19 23 IGNORE WORDS 20 23 INCLUDE 20 20 23 25 53 54 LITERAL 20 48 SECTION 20 20 52 54 SSS 20 52 54 SS 20 52 54 SUBSECTION 20 20 52 54 SUBSUBSECTION 20 20 52 54 S 20 52 54 4 19 54 19 19 19 20 54 21 20 22 31 32 52 53 20 22 53 11 25 absolute permeability 6 33 absTemp 36 absTolIF f 27 addition 21 adsorption 11 advection 9 10 26 28 Qi 6 angled brackets 25 Armijo Goldstein condition 16 29 base unit 21 base units 18 21 26 54 baseLength 26 baseMass 26 baseTemperature 26 baseTime 26 batch system 11 16 bc_f_time 44 bc_f_value 44 bc_t_time 44 bc_t_value 44 bcPresType in 43 bcRegion 43 bcRegion x max 42 bcRegion x min 42 bcRegion y max 43 bcRegion y min 42 bcRegion z max 43 bcRegion z min 43 bcType_f 43 43 bcType t 44 44 boolean 25 bottomhole production well 45 46 boundary condition 43 44 Dirichlet 43 44 inflow 43 44 84 Neumann 43 no flow 43 44 noflow 26 outflow 43 44 Robin 44 BOUNDARY CONDITIONS 42 boundary region 42 branching 40 branching 40 40 branchRatios 40 by 34 34 C 10 40 Ci
12. 3 for information on a data block j For each of the nx 1 ny 1 nz 1 z grid points give the following lt real of a data block L gt zp The z grid points See 84 7 3 for information on a data block 5 5 Porous Medium Material Properties S MATERIAL PROPERTIES real L T 2 9 8 m sec 2 gt gravity Gravitational constant The gravitational constant real grid array 1 gt porosity Porosity A grid array of porosity values given over the 3 D domain as cell centered data one per cell see 4 7 3 An immobile phase can be accounted for by decreasing the porosity of the rock itself and entering those values here instead Herein pore volume refers to flowing phase volume 5 5 1 Permeabilities or Conductivities SS Permeability SS Conductivity This sectional unit is needed only if computeFlow 0 The meaning of the input permeability perm values is either the absolute permeability or the hydraulic conductivity depending on the subsection name The hydraulic conductivity is K pgk po Permeability is assumed internal to the code If an immobile fluid phase is present its effect on the flowing phase should be reflected in the values given below for the permeability perm i e give the product of the absolute permeability and the endpoint relative permeability converted to conductivity if necessary scalar diagonal symmetric face permType scalar diagonal symmetric or
13. 41 lt real grid array L T gt velZ Z velocity A grid array of z Darcy velocity values given over the 3 D domain as point centered data in z and cell centered data in x and y for a total of nx ny nz 1 velocities see 4 7 3 5 12 2 Transport SS Transport This sectional unit is needed only if computeTransport 0 computeChemistry 0 or computeRND 0 For each of the nSpecies species give the following real grid array 1 L 3 conc Molar concentration A grid array of initial molar concentration values given over the 3 D domain as cell centered data one per cell see 4 7 3 The units specification facilities see 4 7 1 can be used to convert mass or other concentrations to molar concentrations 5 13 Boundary Conditions S BOUNDARY CONDITIONS lt integer gt nBCRegions Number of boundary regions The number of distinct regions of the domain boundary lt integer gt maxnBCInterp_f Maximum number of flow interpolation times The maximum number of interpolation times used to specify the flow boundary conditions see below lt integer gt maxnBCInterp_t Maximum number of transport interpolation times The maximum number of interpolation times used to specify the transport boundary condi tions see below The following grid arrays define an order to the boundary regions used elsewhere The following lines appear only if periodicBC_x is 0 lt integer grid array gt bcRegio
14. Currently not used o Il 1 verbosity Driver verbosity flag Boolean to monitor progress in the driver routines A minimal amount of monitoring output will always be given o 111112113114 debug Driver debug flag Flag to set the debug mode in the driver routines The command line argument t gives the options 0 for no debugging tests or output 1 for debugging tests of the parallel set up 2 for debugging output of miscellaneous items 3 for debugging output dump of the input data variables 4 for debugging output dump of the internal global variables The debugging level chosen includes everything at a lower level Level 2 and above produces output Output is written by each processor into its own file called debug where the processor number replaces the question marks word of characters outDir Output directory The directory in which to write output data files Use the LITERAL command if the name begins with a special symbol The current directory is in unix Maximum of 120 characters oll 1 outFormat3D Output 3 D file format Flag selecting the format for 3 D output data files The command line argument t gives the options 1 Eye visualization package format 0 raw data format 1 Tecplot visualization package format lt integer gt nSpeciesPerOutfile Number of species per 3 D output file If multiple species can appear in a single 3 D output file this sets the maximum number allowed oll
15. Number of sorption types The number of distinct sorption or retardation capacities For each of the nSorpTypes sorption types give the following This section induces an ordering of the sorption capacities that is used elsewhere 35 lt real grid array 1 gt porosity sorp Sorption capacities The sorption capacities c The option porosity implies that the retardation is propor tional to porosity Otherwise give one value per cell see 4 7 3 5 6 Phase Properties S PHASE PROPERTIES lt real M LT gt fluidViscosity Flowing phase viscosity The flowing phase viscosity of the resident fluid real M L 3 fluidDensity Flowing phase mass density The average flowing phase mass density The fluid is assumed incompressible so the density of the fluid is not affected by its composition i e we assume dilute solutions This can cause inconsistencies if heavy or light components are present in the fluid in high concentrations integer nPhases Number of phases The total number of fluid and solid phases 5 6 1 Phase Description SS Phase Include one such sectional unit per phase as specified by nPhases This section induces an ordering of the phases that is used elsewhere lt word of characters gt phaseNames Phase name The name of the phase Maximum of 15 characters multi species single species phaseType multi species or single species Whether the phase consists of multiple c
16. Ot oam Qc i by characteristic traking The result is i o At Qc i Then c T 0405 For the Godunov Methods we solve more directly T Ci Cj Mm V uU qei 9 aisi AL Radionuclide decay We solve the equation Oc 6 aos 57 0101 RD C1 sen explicitly by Ci Cj At RP 295 Chemistry We solve for nonequilibrium reactions the equation Oc 6 E RE c m Cn explicitly by i Ci At where At At aici We then allow for equilibrium reactions satisfying an equation of the form RE En RE 61 65 0 Diffusion Dispersion We solve the equation Oc 10 V DVG 6 0401 V DVe 0 implicitly by cU i aioi ADM YAA 9 aici A e 2 3 3 2 Characteristics Mixed Method CMM An explicit characteristics method can be used for advection 2 10 The code scales nearly linearly in parallel 5 No CFL time step constraint is imposed other than that related to the domain decomposition and relatively large time steps can be taken The scheme has minimal numerical dispersion but can suffer from numerical mass and or volume imbalances It is also relatively computationally expensive 2 3 3 3 Higher Order Godunov Method HOG An explicit formally second order Godunov method can be used for advection 15 A postprocess ing step improves the order of accuracy except near sharp fronts or shocks T
17. The Mathematics of Finite Elements and Applications VIII MAFELAP 1993 J R Whiteman ed New York 1994 Wiley pp 199 213 A characteristics mized finite element method for advection dominated transport prob lems SIAM J Numer Anal 32 1995 pp 404 424 T ARBOGAST M F WHEELER AND I YoTOV Logically rectangular mixed methods for groundwater flow and transport on general geometry in Computational Methods in Water Resources X Vol 1 A Peters et al eds Dordrecht The Netherlands 1994 Kluwer Academic Publishers pp 149 156 Logically rectangular mixed methods for flow in irregular heterogeneous domains in Computational Methods in Water Resources XI A A Aldama et al eds vol 1 Southampton 1996 Computational Mechanics Publications pp 621 628 Mixed finite elements for elliptic problems with tensor coefficients as cell centered finite differences SIAM J Numer Anal 34 1997 pp 828 852 82 14 15 16 17 18 19 20 21 22 23 24 25 26 27 30 L C COWSAR C A S SOUCIE AND I YoTOv Parcel v1 04 user guide Tech Rep TICAM 96 28 Texas Institute for Computational and Applied Mathematics The University of Texas at Austin Austin Texas June 1996 C DAWSON Godunov mixed methods for advection diffusion equations in multidimensions SIAM J Numer Anal 30 1993 pp 1315 1332 C DAWSON AND M F WHEELER
18. Transverse dispersion Use this syntax if uniformDispersion is false i e 0 real grid array L 2 T gt molDiffAry Molecular diffusion array lt real grid array L gt longDispAry Longitudinal dispersion array lt real grid array L gt transDispAry Transverse dispersion array SS Sorption This sectional unit is needed only if computeTransport 0 computeChemistry 0 or computeRND 0 lt integer gt nSorpTypes Number of sorption types For each of the nSorpTypes sorption types give the following real grid array 1 gt porosity sorp Sorption capacities S PHASE PROPERTIES 61 real M LT gt fluidViscosity Flowing phase viscosity real M L 3 gt fluidDensity Flowing phase mass density lt integer gt nPhases Number of phases SS Phase Include one such sectional unit per phase as specified by nPhases lt word of characters gt phaseNames Phase name multi species single species phaseType multi species or single species real 1 L 3 lt omit gt chmPhaseDensity Molar density S CHEMISTRY This sectional unit is needed only if computeChemistry gt 0 real temperature reading gt absTemp Equilibrium temperature real M LT 2 gt chemPres Equilibrium pressure ideal non ideal ideal ideal or non ideal S CHEMICAL SPECIES integer nComps Number of components integer nProducts Number of products nSpecies
19. a larger number of Lagrange multi pliers The UNSF is prohibitively slow for realistic simulations and is included mainly for archival reasons The choice of version is made at the time of compilation 2 4 3 Numerical Solution An interior point optimization technique 19 26 is used to solve for the minimum of the Gibbs free energy 2 5 General Chemistry We reproduce here for convenience information from the chemistry manual 27 Basic reaction theory is presented The variables in typewriter font refer to variables to be input as data as described in the chapter on the input data file s Chapter 5 2 5 1 Thermodynamics In this section we discuss the issues of phase identity of a species the usage of practical concen tration scales versus computationally convenient ones and unit conversions The versions that are based on mole fractions i e SF and UNSF differ from the molar concentration based version RNSF and we discuss these cases separately In this section we use notation established by Saaf in his Ph D Thesis 25 In particular the system is comprised of Ns species nSpecies of which Nc are components nComps and Ng products nProducts When convenient we distinguish N kinetic and Ng equilibrium products Ng N ES N Finally we let Nm denote the number of minerals and Im the corresponding index set 11 2 5 1 1 Molar concentration based system version RNSF When the RNSF is used the composition of the
20. absence of minerals If minerals are present this set is augmented by one complementarity and one optimality equation for each mineral It can be shown that mineral concentrations in the computed solution will satisfy the inequality No chmAbsTol 1 c7 log K c 5 chmAbsTol ln 23 Neti J j l which is a satisfying result In particular if the mineral is present at a large concentration cy 4 7 1 this is approximately Nc K ej 1 chmabsTo1 2 22 j l For numerical reasons a non zero tolerance chmEpsConc gt 0 is used in enforcing the condition of non negativity of species When the SF is used all species concentrations affected by the equili bration are bounded below by chmEpsConc For the RNSF which uses logarithmic transformations of some of the variables this restriction is only in place for mineral species other species can be arbitrarily close to zero It is assumed that chmEpsConc lt chmAbsTol The variable chmTestSolFlag enables a test of the second order sufficiency conditions at the computed solution i e that the Hessian projected onto the tangent space of the active constraints have all positive eigenvalues see e g 30 This procedure can only be used with versions UNSF and SF It is probably of limited value for most simple aqueous ideal systems which have known convexity properties It could prove useful for batch testing of more complicated multi phase systems 2 5 3 2
21. cell of region 2 2 4 histIndexHi_t_x _y _z Upper cell of region 2 nHistSpecies_t Number of species in the history 1 2 32 histSpecie_t Species in the history dStep0ut_hist_t Steps between history outputs 1 dtOut_hist_t Time between history outputs A 3 Sample input file 2 A 3 1 The main input file SAMPLE GENERAL INPUT FILE FOR PARSSIM1 version 2 1 S GENERAL INFO 1 computFlow 2 computeTransport 1 computeChemistry O computeRND m baseLength Internal length base units g baseMass Internal mass base units hr baseTime Internal time base units degC baseTemperature degK degC or degF S SOLUTION PARAMETERS 1 1 1 nDom x y _z Parallel subdomain divisions 1 auto select SS Flow 100 maxIterIF f Interface maximum number of iterations 1 0e 6 relTolIF f Interface relative tolerance 1 0e 12 absTolIF f Interface absolute tolerance 5 pcTypeIF_f Interface preconditioner 2048000 dWorkspace_f Double workspace for flow SS Transport 100 maxlIter t Dispersion maximum number of iterations 1 0e 6 relTol t Dispersion relative tolerance 204800 dWorkspace t Double workspace for transport 74 1 cflFactor Numerical CFL factor SS Chemistry 0 chmLoadBalFlag Parallel chemistry load balancing 100 chmMaxIter Maximum number of nonlinear iterations 1 0E 09 chmAbsTol Absolute KKT tolerance 1 0e 30 chmEpsConc Minimal concentrati
22. constant 22 22 coordinate direction 6 31 copyright 7 D 6 Darcy velocity 6 7 41 42 data block 22 23 31 32 53 data items 19 21 debug 48 52 debug 48 degC 21 26 degF 21 26 degK 21 26 At 6 density 6 depth 31 diagonal 33 diagonal 33 34 diffusion dispersion 9 10 10 26 27 50 diffusion dispersion tensor 6 8 35 dilute solutions 7 12 13 36 Dirichlet condition 43 44 Dispersion 35 dispersion see diffusion dispersion dispersion tensor see diffusion dispersion tensor dispersive flux 44 dissolution 11 distance from equilibrium 11 14 15 division 21 domain decomposition 5 7 10 49 Draw 50 driver 5 81 dStep ut conc 49 dStep Out hist t 51 dStepOut pres 49 dStep ut vel 49 dtMax f 30 dtMax t 28 31 dtOut_conc 49 dtOut_hist_t 51 dtOut_pres 49 dtOut_vel 49 dWorkspace_c 29 dWorkspace_f 27 dWorkspace t 27 echo 18 echoChem 18 end of line 19 epsConc 17 equilibrium 5 38 error messages 52 exponentiation 21 extraction well 8 45 Eye 5 48 50 56 57 face 33 face 33 35 face centered data 22 35 57 Fahrenheit 21 First Order Godunov Method 10 flag 25 Flow 27 30 41 43 46 49 flow 5 5 6 20 81 flow equation 6 fluidDensity 6 36 fluidViscosity 6 7 36 39 40 FOG see Godunov Method Fortran 10 22 27 29 30 40 57 function of time 43 43 44 44 46 47 g 6 GENERAL INFO 25
23. cuts The number of times the transport time step size can be cut per step attempted integer 2 integer 2 integer 2 gt padx pady padz Number of cells for subdomain overlap pad The number of extra grid cells in each coordinate direction in the overlap region padded or shadow region on each end of each subdomain lt real 1 10 gt volTol Volume discrepancy tolerance The tolerance for the transport s relative volume discrepancy lt integer 1 gt lt integer 1 gt lt integer 1 gt volRefine_t_x volRefine_t_y volRefine_t_z Trace back volume refinement Trace back volume refinement in each direction The value 1 is no refinement full grid cells will be traced back otherwise each grid cell can be divided along the coordinate directions into an integral number of pieces and each sub cell will be traced 5 2 3 Chemistry SS Chemistry This sectional unit is needed only if computeChemistry gt 0 O l 1 1 gt chmLoadBalFlag Parallel chemistry load balancing Specifies whether automatic parallel load balancing of the computational work among proces sors should be attempted This creates parallel communication overhead but generally speeds up the computation in parallel lt integer 100 gt chmMaxIter Maximum number of nonlinear iterations The maximum number of iterations allowed in the nonlinear solver lt real 1 gt chmAbsTol Absolute KKT tolerance The absolute tolerance on the residual
24. enables one to include the numerical values of unsupported units It is relatively easy for a programmer to add a new unit name to the code or to change the initial base units used internally in the code Temperature degree increments are an ordinary case of units as described above however temperature readings are a special case These need simply the temperature scale used degrees Kelvin degK Celsius degC or Fahrenheit degF Examples follow 21 2e 5 cm sec 2 floating point notation is allowed 8 2 kg m sec 2 negation and integral exponentiation are allowed 4 37 2 dimensions may be incorrect value is doubled 50 yr initial division is allowed 1 2 2 1e 4 6 8 1 cm m 2 3 cm sec 4 nonsense but a correctly formed expression dk dk HOH OH OF 100 degC a temperature reading 4 7 2 Data Blocks and the Repetition Symbol Data items may come in a block of data as when a grid array see below is specified Such a data block is bounded by the list delimiter symbols and as described above Within such a block the repetition character may be used A set of n single data items each with the value value may be indicated by n value Proper white space may precede value but only spaces may precede e Units may be specified for the real numbers in a data block but only for the list as a whole that is only a single expression can be given and
25. face Declare the permeability to be a scalar one value per cell diagonal tensor three values per cell a full symmetric tensor six values per cell or face centered diagonal tensor values a single value for each face of the grid This must be followed by one of the following sets of input data items 33 5 5 1 1 Scalar Permeabilities Use the following when permType is scalar lt real grid array L 2 permeability or L T conductivity gt perm Scalar permeabilities A grid array of scalar permeability values given over the 3 D domain as cell centered data one per cell see 4 7 3 5 5 1 2 Diagonal Tensor Permeabilities Use the following when permType is diagonal by cells by components permGrouping Group by grid cells or by tensor components Give the entire permeability tensor for each grid cell or give successively a single component of the tensor for the entire grid The word is optional xx yy zz lt any permutation gt permComponentOrder Order of tensor s components Declare the order of the permeability tensor component data lt 1 or 3 real grid arrays L 2 permeability or L T conductivity gt perm Diagonal tensor permeabilities One or three grid arrays of diagonal permeability values given over the 3 D domain as cell centered data three or one number per cell see 4 7 3 If the permGrouping is cells give the three components of the tensor as ordered by permComp
26. gt padx pady padz 58 Number of cells for subdomain overlap pad lt real 1 10 gt volTol Volume discrepancy tolerance lt integer 1 gt lt integer 1 gt lt integer 1 gt volRefine_t_x volRefine_t_y volRefine_t_z Trace back volume refinement SS Chemistry This sectional unit is needed only if computeChemistry gt 0 O 1 lt 1 gt chmLoadBalFlag Parallel chemistry load balancing lt integer 100 chmMaxIter Maximum number of nonlinear iterations real 1 gt chmAbsTol Absolute KKT tolerance real moles L 3 1 0e 16 gt chmEpsConc Minimal concentration parameter O 1 lt 0 gt chmScaleFlag Diagonal scaling O 1 lt 0 gt chmTestSolFlag Test second order sufficiency conditions 1 0 1 lt 1 gt chmGuessType initial guess 1 auto 0 transported 1 previous O 1 lt 0 chmInterpFlag Use interpolation in the back track line search lt integer 10 gt nViol Number of non monotonic line searches real 1 1 0e 4 gt chmAlpha Alpha parameter O 1 1 gt hessFlag Analytical Hessian real 1 0 8 gt tauMin Movement to the boundary factor real 1 1 0e 2 gt chmRho IP reduction factor of perturbation parameter integer ntReact Number of reaction steps per transport step O 11 1 I 2 odeAlgType ODE integration 0 RK1 1 RK2 2 RK4 O 1 lt 0 gt switchFlag Species switching lt integer 1000 dWorkspace_c Double
27. it applies to each number of the data block The units expression must appear just before the list of numbers as in cm 3 sec 2 1 3 2 407 5 4 7 3 Grid Arrays Constant Nearly Constant and Fully Specified Data may be distributed in space and therefore assigned to every cell or point of the grid in 3 D or in 2 D to every face or point of the domain boundary The data per grid cell may be a single numerical value or a set of numerical values of a given fixed size A uniform syntax is used to specify such an array Within the array the x coordinate varies most rapidly then y and finally z as in Fortran For cell centered data the index range goes from 1 to the number of grid cells in the given direction i e the full range is 1 nx 1 ny 1 nz For point centered data the index range goes from 0 to the number of grid cells in the given direction i e the full range is 0 nx 0 ny 0 nz For face centered data the index range depends on the type of face considered For an x face the array is point centered in the x direction and cell centered in the other directions similarly for a y face and a z face i e the full range is 0 nx 1 ny 1 nz 1 nx 0 ny 1 nz and 1 nx 1 ny 0 nz for x face y face and z face face centered data respectively A grid array that consists of a single set of repeated values i e a constant set of values can be specified as constant lt value s gt where if the values are real number
28. nComps nProducts SS Component SS H Component SS H20 Component SS Product Include one such sectional unit per component as specified by nComps Follow by one such sectional unit per product as specified by nProducts lt word of characters gt specieName Specie name real M molecularWeight Molecular weight lt real 1 gt phaseDist The alpha parameter in Henry s Law integer lt omit gt sorpType Sorption type SSS Chemistry This sectional unit is needed only if computeChemistry gt 0 Use this syntax if the subsection heading is Component H Component or H20 Component lt integer gt phaseldentity Participating phase real 1 gt J compCharge Charge Use this syntax if the subsection heading is Product lt integer gt phaseldentity Participating phase For each of the nComps components give the following real 1 stoich Stoichiometry formula number O II 14 I 2 reactionType Reaction type Use this syntax if reactionType specifies an equilibrium reaction real 1 gt pK Log 10 equilibrium constant Use this syntax if reactionType specifies a mass action type kinetic reaction real 1 gt real 1 gt pKf pKb Log 10 forward and backward rate constants For each of the nComps components give the following lt real 1 gt stoichK Rate law powers J 62 Use this syntax if reactionType specifies a monod type kinetic reacti
29. number of iterations Maximum number of iterations for solving the transport diffusion dispersion problem lt real 1 1 0e 6 gt relTol_t Dispersion relative tolerance Relative tolerance for solving the transport diffusion dispersion problem lt integer 10000 gt dWorkspace_t Double workspace for transport The amount of Fortran double precision workspace to allocate for the transport routines The transport routines will abort if this value is too small In that case raise its value and rerun the code 27 lt real 1 1 gt cflFactor Numerical CFL factor The numerical CFL constraining factor applies to the advection schemes If this number is negative the desired fixed time step size dtMax_t will be attempted If this number is positive an automatic selection will be made If the CMM option is selected by computeTransport then the CFL limit refers to the maximum amount of characteristic tracking that may take place in each subdomain s overlap region padded or shadow region specified by padx pady and padz The cflFactor will then be used to select a time step that is the given fraction cflFactor of the maximum estimated time step that would keep each subdomain characteristic track within its subdomain or overlap region 5 2 2 1 Characteristics Mixed Method SSS_ CMM This sectional unit is needed only if computeTransport selects the CMM advection scheme integer 2 gt ntCutMax t Maximum number of time step
30. ny integers integer grid array bcRegionz max Boundary region map A grid array of boundary region identification numbers from 1 to nBCRegions given over the 2 D boundary face for maximal z This is cell centered data one per cell face see 84 7 3 for a total of nx ny integers 5 13 1 Boundary Region Specification SS Region Include one such sectional unit per boundary region as specified by nBCRegions These regions are mapped to the physical boundary and ordered by the bcRegion grid arrays 5 13 1 1 Flow SSS Flow This sectional unit is needed only if computeFlow 0 oll itl 2 bcType f BC type Flag declaring the type of the boundary condition for flow The command line argument t gives the options 0 for no flow condition u v 0 0 normal velocity 1 for a Neumann condition u v qn specified normal velocity where the inflow condition is positive and the outflow condition is negative 2 for a Dirichlet condition p py specified pressure Above gn or py needs to be specified as a function of time Function of Time Specification The following lines appear only if a function of time needs to be specified that is bcType f 0 lt omit gt pressure potential head hydraulicHead bcPresType_in Dirichlet pressure potential head or hydraulicHead If bcType f selects the Dirichlet condition give the meaning of the variable to which the condition applies see page 49 for a
31. or chemistry without modifying the input file s If there is any doubt one can use the command line argument e or E to see if a given sectional unit was actually processed 4 5 List Delimiters List delimiters enclose lists that have a user specified number of items A list is begun with a begin list symbol and ended with an end list symbol P The delimiters are used to check that the user has counted the number of list items correctly 20 4 6 Key words Key words are specified strings of characters Spaces may be embedded in a key word 4 7 Data Items A data item can be an integer not be followed by a period a real number in standard floating point or exponential notation a word of characters terminated by a space newline or tab or a line of characters including spaces and tabs terminated by a newline character Data items of characters are limited in length The maximum allowed number of characters is 120 for a directory name 50 for a file name 120 for a line of character 30 for a word of characters 15 for phase and species names 4 7 1 Physical Units Real numbers have physical units Physical units can be optionally assigned by using the square bracket notation lt real gt lt units expression gt where the units expression may contain integer and real numbers the names of units arithmetic operators parentheses and spaces which are ignored The statement declares the physical uni
32. per well as specified by nWells O DT X bp 2 ES Pod wellType Well type Flag specifying the type of the well T he command line argument t gives the options 0 for an inactive well 1 for an injection well 2 for an extraction well or production well 3 for an bottomhole production well that uses the Peaceman correction 23 4 for an reinjection well or sparging well A reinjection well has both an extraction and an injection interval the fluid extracted is reinjected elsewhere An inactive well is equivalent to no well at all and no further specification should be given below If the well is not inactive include the following information real L gt wellRadius Well radius The radius of the well integer integer hWellIndex The x and y indices of well The location of the well in the grid given as the x and y cell centered data index numbers The well is generally assumed to occupy entire cells integer zWellIndexLo Lower z index of well The top if zIsDepth is true or bottom otherwise of the well in the grid given as the lower z cell centered data index number If wellType selects a reinjection well specify the injection part integer zWellIndexHi Higher z index of well The bottom if zIsDepth is true or top otherwise of the well in the grid given as the upper z cell centered data index number If wellType selects a reinjection well specify the injection part omit integer
33. reactions respectively see 2 5 3 3 oO 1 lt 0 switchFlag Species switching Boolean to enable the species switching procedures see 2 5 3 3 lt integer 1000 gt dWorkspace_c Double workspace for chemistry The amount of Fortran double precision workspace to allocate for the chemistry routines The chemistry routines will abort if this value is too small In that case raise its value and rerun the code The amount of workspace needed depends on the number of species and the reactions but not on the grid size 29 lt integer 1000 gt iWorkspace_c Integer workspace for chemistry The amount of Fortran integer workspace to allocate for the chemistry routines The chemistry routines will abort if this value is too small In that case raise its value and rerun the code The amount of workspace needed depends on the number of species and the reactions but not on the grid size 5 3 Time Control S TIME real T 0 tInitial Initial time The initial or starting time real T gt tFinal Final time The final or ending time 5 3 1 Flow SS Flow This sectional unit is needed only if computeFlow 0 integer 4 nStepsMax f Maximum number of flow steps The maximum number of flow steps allowed integer ndtMax f Number of flow time step sizes The maximum number of flow step size increments given below For each of the ndtMax_f time step sizes give the following The data are the
34. system is computed internally in molar concentra tions cj conc The standard mass action expressions apply that is Nc enoti TT i 1 Ng 2 1 j 1 with a special case for minerals i Im Nc Kill Gy X ifcenoti gt 0 2 2 j l No K I cy lt 1 if CNoti 0 j l where K is the equilibrium constant for the ith reaction and aj are stoichiometric coefficients stoich In Parssim1 all species participating in multi species phases including species in the aqueous phase sorbed phases or ion exchange sites should have phaseIdentity 1 since the chemical potential model is the same for these species based on molar concentrations Minerals however constitute new phases and should be labeled with phaseIdentity 2 Nm 1 The use of the water component is optional in the RNSF If water is not included it is an implicit component see e g 21 If it is included an additional mass balance equation for water will result Note that the molar phase density chmPhaseDensity should be set to a representative value for the aqueous phase approximately 55 4 mol liter if water is included as a component 2 5 1 2 Mole fraction based system versions SF and UNSF The versions SF and UNSF are based on internally computing chemical equilibrium in terms of mole fractions This treatment is more general than the use of molar variables but it also carries some limitations In particular this formulation allows
35. the format used by the visualization tool Eye or Tecplot or raw data can be written Chapter 2 Governing Equations and Numerical Methods We describe briefly the flow transport chemistry and radionuclide decay portions of the simulation Each is solved independently of the others by making use of time splitting techniques 16 17 as described later in this chapter Some general references to our approach can be found in 9 29 and in a multi phase setting in 7 1 6 See also the general reference 22 2 1 Notation We use the following standard symbolic notation throughout this manual cj is the molar concentration of species given in moles per pore volume i e flowing phase volume conc D isthe diffusion dispersion tensor g is the gravitational constant gravity k is the absolute permeability perm p isthe pressure q isa source or sink i e it represents the wells and leaking sources t is time u is the Darcy velocity velX velY ve1Z wi is the molecular weight of species nolecularWeight x is the first logical coordinate direction y is the second logical coordinate direction z isthe third logical coordinate direction a is the sorption coefficient essentially the Henry s Law constant phaseDist At isthe current time step size u is the flowing phase viscosity Ho is the viscosity of the resident fluid fluidViscosity v isthe outward unit normal vector to the domain p is the flowing pha
36. the simulation of multi species phases other than the aqueous phase such as a gas phase This can be accomplished by listing in additional to the aqueous phase other multi species phases in the input file s As before minerals are treated as separate phases and should correspond to different values of phaseIdentity A difficulty arises in using practical concentration units with a mole fraction based calculation To be compatible with the use of molar concentrations in aqueous chemistry a conversion of the equilibrium data is performed in the code The conversion modifies the equilibrium constants based on the assumption of dilute solutions and the specified value of chmPhaseDensity for each multi species phase We illustrate this conversion for the simple case of one aqueous multi species phase The molar concentration of a species c is defined as n n nea 2 Cl V a V 3 where V is the aqueous phase volume and the total number of moles in the aqueous phase Introducing the molar phase density peoa I 2 4 p V 2 4 a and using the definition of the mole fraction of the lth species we can rewrite a RE 2 5 Substituting this in the original mass action expression 2 1 yields Nc x ji ym INoti Kj s peg 2 6 j l where the modified equilibrium constant A7 satisfies log K e aji log p 4 if i is a mineral 2 7 log K xm aji l logp otherwise vum log K7
37. time it begins and its size The first time step size will be assumed to be valid at the initial time The times must increase monotonically lt real T gt tdtMax f Time of onset Time of onset of the next time step size real T gt dtMax f Time step size The maximum allowable time step size 5 3 2 Transport SS Transport This sectional unit is needed only if computeTransport 0 computeChemistry 0 or computeRND 0 lt integer gt nStepsMax t Maximum number of transport steps The maximum number of transport steps allowed lt integer gt nStepsMax t per f Maximum number of transport steps per flow step The maximum number of transport steps allowed per flow step A negative number is in trepreted as infinity integer ndtMax t Number of transport time step sizes The maximum number of transport step size increments given below 30 For each of the ndtMax_t time step sizes give the following The data are the time it begins and its size The first time step size will be assumed to be valid at the initial time The times must increase monotonically lt real T gt tdtMax t Time of onset Time of onset of the next time step size lt real T gt dtMax t Time step size The maximum allowable time step size 5 4 Spatial Grid S GRID lt integer gt lt integer gt lt integer gt nx ny nz Number of grid cells The number of grid cells i e elements in each coordinate direction
38. used within the code Note however that the base units can be modified by the user v Just display Parssim1 version information h Print a summary of the above usage information only 18 Chapter 4 General Input File Syntax and Physical Units All input is read from one or more data files Key words and data items are arranged in sections subsections and sub subsections and ordered item by item possibly separated by white space comments and list delimiters and possibly special commands 4 1 White space White space is used to separate other items and is otherwise ignored on input White space consists of the characters for space tab newline or end of line comma semicolon and colon White space is ignored on input The punctuation characters allow the user to delineate sets of items for example Proper white space consists of the characters for space tab and newline Command names and 5 key words are terminated by proper white space 4 2 Commands Commands alter the way the data file s are read Commands take the form command symbol followed immediately no space by a command name Note that the command symbol is ignored in any explicit comment and in skipped sectional units but not before the first section see the discussion below on comments and sectinal units A list of commands follows Null command Null command Equivalent to BEGIN COMMENT 4 Equivalent to
39. with this choice the reaction rate is zero if and only if the reaction satisfies the equilibrium conditions with an equivalent equilibrium constant of the form k 2 14 This follows from the fact that K 0 lt gt kf l ri ays m Kenoti 0 j 1 E lt 2 c No i 0 i j l No A lt CNo i K TOKE j 1 The definition of the equilibrium constant for the completed reaction 2 14 allows species switching switchFlag to be used see also 2 5 3 3 We point out that regardless of the choice of equilibrium module the rate laws always use molar concentrations 2 5 2 2 Monod style expressions Monod introduced an empirical expression of the form dX S e oue ee dt K S to account for the growth of microbes X on a substrate S The constant K is commonly referred to as the half saturation constant halfSatConst This expression is widely used to model biologically mediated reactions We allow a generalized version of this expression 2 15 No No n k e T Sep nets Gee NE 2 16 j 1 j1 Ky sp 14 where K ae are the half saturation constants for the component species more than one component is allowed to participate in each reaction and N M is the number of Monod style rate expressions Any or all or none of the N ri kinetic reactions may have Monod style expressions i e 0 lt N a lt N K When RU 0 for all j the Monod expression reduces to the distance from equilibrium express
40. workspace for chemistry lt integer 1000 iWorkspace c Integer workspace for chemistry S TIME real T 0 tInitial Initial time real T gt tFinal Final time SS Flow This sectional unit is needed only if computeFlow 0 integer nStepsMax f Maximum number of flow steps integer ndtMax f Number of flow time step sizes 1 For each of the ndtMax f time step sizes give the following real T tdtMax f Time of onset real T gt dtMax f Time step size SS Transport This sectional unit is needed only if computeTransport 0 computeChemistry 0 or computeRND 0 lt integer gt nStepsMax_t Maximum number of transport steps integer nStepsMax_t_per_f Maximum number of transport steps per flow step lt integer gt ndtMax_t Number of transport time step sizes For each of the ndtMax_t time step sizes give the following lt real T gt tdtMax_t Time of onset lt real T gt dtMax_t Time step size S GRID lt integer gt lt integer gt lt integer gt nx ny nz Number of grid cells 59 lt integer gt lt integer gt lt integer gt periodicBC_x _y _z Boolean for periodicity O 1 zIsDepth Direction of the z coordinate uniform rectangular xy rectangular 4 gridType uniform rectangular xy rectangular or Use the following when gridType is uniform real L gt real L gt xMin xMax Minimal and maximal x p
41. 0 min tFinal Final time SS Flow 1 nStepsMax f Maximum number of flow steps 1 ndtMax f Number of flow time step sizes 0 tdtMax f Time of onset 1 hr dtMax f Time step size SS Transport 2 nStepsMax_t Maximum number of transport steps 1 nStepsMax t per f Maximum number of transport steps per flow step 1 ndtMax t Number of transport time step sizes 1 0 tdtMax t Time of onset 30 s dtMax t Time step size 8 GRID 8 4 4 nx ny nz Number of grid cells 0 0 0 periodicBC x Jy _z Boolean for periodicity 1 zIsDepth Direction of the z coordinate uniform gridType uniform rectangular xy rectangular or 0 18 m xMin xMax Minimal and maximal x points 68 O 12 n yMin yMax Minimal and maximal y points 1 100 cm zMin zMax Minimal and maximal z points 8 MATERIAL PROPERTIES 9 8 m sec 2 gravity Gravitational constant 6400 25 6400 5 porosity Porosity SS Permeability diagonal permType scalar diagonal or symmetric by components permGrouping Group by grid cells or by tensor components XX yy ZZ permComponentOrder Order of tensor s components perm Diagonal tensor permeabilities md 40 45 47 43 40 45 47 43 39 36 2038 39 36 2038 39 43 41 42 39 43 41 42 45 47 55 65 45 47 55 65 50 55 57 53 50 55 57 53 49 46 2048 49 46 2048 49 53 51 52 49 53 51 52 55 57 65 75 55 57 65 75 40 45 47 43 40 45 47 43 39 36 2038 39 36 203
42. 1 initial0ut Output 3 D initial conditions flag Boolean for giving output data files at the initial time 1 outFormatiD Output 1 D file format Flag selecting the format for 1 D output data files The command line argument t gives the options 1 Tecplot visualization package format 48 5 16 1 Flow SS Flow This sectional unit is needed only if computeFlow 0 lt integer gt dStep ut pres Steps between pressure outputs The number of flow steps between writing pressure output data files The value 1 is interpreted as positive infinity i e generate no output due reaching a certain step number integer dStep ut vel Steps between velocity outputs The number of flow steps between writing velocity output data files The value 1 is interpreted as positive infinity i e generate no output due reaching a certain step number real T gt dtOut_pres Time between pressure outputs The flow time increment between writing pressure output data files The value 1 is interpreted as positive infinity i e generate no output due reaching a certain time real T gt dtOut vel Time between velocity outputs The flow time increment between writing velocity output data files The value 1 is interpreted as positive infinity i e generate no output due reaching a certain time potential pressure hydraulicHead head presType out potential pressure hydraulicHead or head The physical interpreta
43. 1111111111111111111100 0011111111111111111111111111111111111100 0011111111111111111111111111111111111100 Q O111114117711021410110 11114111127212427y412114731141111 900 0011111111111111111111111111111111111100 Q 0 1 1 111111 111 141411112111 1111111 1111411111 900 0011111111111111111111111111111111111100 0011111111111111111111111111111111111100 0011111111111111111111111111111111111100 0011111111111111111111111111111111111100 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 80 Authors and Acknowledgments Parssiml was developed in the Center for Subsurface Modeling Many people contributed to its development The primary authors are acknowledged below Driver The driver and overall code framework was developed primarily by Todd Arbogast Contribu tions were made also by Ashokkumar Chilakapati Lawrence C Cowsar Robert McLay Douglas Moore and Ivan Yotov Flow Parcel 14 was developed initially by Lawrence Cowsar and modified by Fredrik Saaf Carol San Soucie Ivan Yotov and Robert McLay Frederic d Hennezel and Todd Arbogast developed the calling and set up routines fo
44. 6 cell 34 cell centered data 22 33 34 41 43 45 47 50 56 57 cell centered finite difference 5 7 10 cells 34 cells 34 Celsius 21 Center for Subsurface Modeling 5 81 CFL 10 28 cflFactor 28 28 Chain 40 chainLen 40 40 chainType 40 40 Characteristics Mixed Method 5 9 10 26 28 50 CHEMICAL SPECIES 37 Chemistry 26 36 37 chemistry 5 5 6 8 11 11 17 20 26 28 36 37 50 81 chemPres 36 chmAbsTol 15 16 28 28 chmAlpha 16 29 chmEpsConc 16 17 28 chmGuessType 16 29 chmInterpFlag 16 29 chmLoadBalFlag 15 26 chmMaxIter 15 28 chmPhaseDensity 12 13 36 chmRho 17 29 chmScaleFlag 29 chmTestSolFlag 16 29 CMM see Characteristics Mixed Method CMM 28 colon 19 25 comma 19 command 19 19 54 command line argument 18 18 20 21 25 27 29 37 39 43 45 48 52 54 55 comment 18 19 20 25 compCharge 38 complexation 11 Component 37 component 34 37 38 components 34 computational domain 5 computeChemistry 26 28 30 35 37 39 40 42 49 computeFlow 25 26 27 30 33 39 41 43 46 49 computeRND 26 30 35 39 40 42 49 computeTransport 26 27 28 30 35 39 42 44 46 49 50 conc 6 8 12 42 50 conc 56 concentration 6 11 42 47 concHist 57 concHist CP 57 concHist CT 57 concHist MP 57 concHist M T 57 concTime 56 Conductivity 33 conductivity 33 conservation of mass 6
45. 8 39 43 41 42 39 43 41 42 45 47 55 65 45 47 55 65 50 55 57 53 50 55 57 53 49 46 2048 49 46 2048 49 53 51 52 49 53 51 52 55 57 65 75 55 57 65 75 constant 40 md nearly constant 10 md 1 1 4 1 4 1 1 20 md SS Dispersion 1 uniformDispersion Flag for nonuniform dispersion 1 cm 2 day molDiff Molecular diffusion 10 cm longDisp Longitudinal dispersion 1 cm transDisp Transverse dispersion SS Sorption 2 nSorpTypes Number of sorption types porosity sorp Sorption capacities constant 50 S PHASE PROPERTIES 2 cp fluidViscosity Flowing phase viscosity O 9 g cm 3 fluidDensity Flowing phase mass density 2 nPhases Number of phases SS Phase 1 OIL_PHASE phaseNames Phase name 69 multi species phaseType multi species or single species 1 cc chmPhaseDensity Molar density SS Phase 2 ROCK phaseNames Phase name Single species phaseType multi species or single species S CHEMISTRY 120 degF absTemp Equilibrium temperature 3 5 atm chemPres Equilibrium pressure ideal ideal ideal or non ideal S CHEMICAL SPECIES 2 nComps Number of components 1 nProducts Number of products SS H Component 1 H specieName Specie name 1 g molecularWeight Molecular weight E phaseDist The alpha parameter in Henry s Law 1 sorpType Sorption type SSS Chemistry 1 phaseIdentity Participating phase 1 compCharge Charge SS Componen
46. Algorithmic parameters for the nonlinear iteration The parameter chmGuessType specifies the type of initial guess used to initialize the nonlinear equilibration iteration The group of parameters chmInterpFlag nViol and chmAlpha relate to the line search glob alization of the nonlinear solver The merit function is the square of the Lz norm of the perturbed KKT conditions see 18 25 26 for details The well known Armijo Goldstein condition see 18 is imposed at each step with the setting nViol 0 if a value greater than zero is specified nViol violations of this condition are accepted before it is again required to hold The parameter chmAlpha is the fraction of the initial decrease that the decrease resulting from the adjusted step is required to have and chmInterpFlag indi cates whether a simple reduction of the step length or a more costly cubic quadratic interpolation procedure is used in the line search step see 18 Often the algorithm converges significantly faster if a few violations of the Armijo Goldstein condition are tolerated and for this reason we recommend values of nViol in the range of 1 to 5 However for highly sensitive problems it may be necessary to use a value of 0 to avoid jeopardizing global convergence It is possible and sometimes useful to combine the setting nViol 0 with chmAlpha 0 This is another way of relaxing the Armijo Goldstein condition by allowing a slight increase in the merit functi
47. An operator splitting method for advection diffusion reaction problems in MAFELAP Proceedings VI J A Whiteman ed Academic Press 1988 pp 463 482 T ime splitting methods for advection diffusion reaction equations arising in contaminant transport in Proceedings ICIAM 91 R O Malley ed SIAM 1992 pp 71 82 J DENNIS AND R B SCHNABEL Numerical Methods for Unconstrained Optimization and Nonlinear Equations Prentice Hall 1983 A S EL BAKRY R A Tapia T TSUCHIYA AND Y ZHANG On the formulation and the ory of the Newton interior point method for nonlinear programming Journal of Optimization Theory and Applications 89 1996 pp 507 541 R GLOWINSKI AND M F WHEELER Domain decomposition and mized finite element meth ods for elliptic problems in First International Symposium on Domain Decomposition Methods for Partial Differential Equations R Glowinski et al eds SIAM Philadelphia 1988 pp 144 172 F M MOREL AND J G HERING Principles and Applications of Aquatic Chemistry Wiley 1993 J C PARKER Multiphase flow and transport in porous media Reviews of Geophysics 27 1989 pp 311 328 D W PEACEMAN Interpretation of well block pressures in numerical reservoir simulation Society of Petroleum Engineers Journal 1978 pp 183 194 T F RUSSELL AND M F WHEELER Finite element and finite difference methods for con tinuous flows in porous media in The Mathematics of Reservoir Simul
48. Armijo Goldstein con dition of monotonicity see 18 is imposed on the step If positive we tolerate that many consecutive violations of this condition before it is imposed again see 2 5 3 2 lt real 1 1 0e 4 gt chmAlpha Alpha parameter The alpha parameter in the Armijo Goldstein condition see 82 5 3 2 and 18 O 1 1 gt hessFlag Analytical Hessian Boolean to select use of the analytic Hessian Otherwise a finite difference computation of the Hessian is performed for the SF and UNSF versions Currently the finite difference code is redundant but included for future use with non ideal systems It is best to use the analytic Hessian option 1 lt real 1 0 8 gt tauMin Movement to the boundary factor Smallest allowable percentage of movement to the boundary in the interior point algorithm see 2 5 3 2 and 25 lt real 1 1 0e 2 gt chmRho IP reduction factor of perturbation parameter The multiplicative reduction of the interior point perturbation parameter see 2 5 3 2 Must be less than 1 lt integer gt ntReact Number of reaction steps per transport step Target number of reaction sub time steps per transport time step The value chosen depends on the relative time scales involved in transport and reactions see 2 5 3 3 oll 1 1 2 odeAlgType ODE integration O RK1 1 RK2 2 RK4 Flag to select first 0 second 1 or fourth order 2 in time accurate Runge Kutta integration of kinetic
49. Brief Summary Parssim1 is a code developed by the Center for Subsurface Modeling CSM of the Texas Institute for Computational and Applied Mathematics TICAM at The University of Texas at Austin Austin Texas It is an aquifer or reservoir simulator for the incompressible single phase flow and reactive transport of subsurface fluids through a heterogeneous porous medium of somewhat irregular geometry It is also capable of simulating the decay of radioactive tracers or contaminants in the subsurface linear adsorption wells and leaking sources and bioremediation Although the code uses very simple rectangular data structures for efficiency and accuracy the subsurface domain can be of irregular geometry The subsurface domain is assumed to be described by a logically rectangular grid A mapping technique is used to map the irregular physical domain with its hills valleys internal faults and strata etc to the rectangular computational domain without loss of accuracy or efficiency The grid may cover a domain that is periodic in one or more coordinate directions i e periodic boundary conditions are allowed The code can run in serial on a single processor or on a massively parallel distributed memory computer or collection of computers using MPI computational results indicate that it has very good parallel scaling properties We use domain decomposition to compute in parallel The grid is divided into subdomains one for each paralle
50. D output file O 1 initialOut Output 3 D initial conditions flag 1 outFormatiD Output 1 D file format SS Flow This sectional unit is needed only if computeFlow 0 lt integer gt dStep ut pres Steps between pressure outputs lt integer gt dStep ut vel Steps between velocity outputs lt real T gt dtOut_pres Time between pressure outputs real T dt ut vel Time between velocity outputs potential pressure hydraulicHead head presType out potential pressure hydraulicHead or head integer 0 1000 gt outFlags f 1 Level of flow verbosity integer 0 1000 outFlags f 2 Level of flow debugging SS Transport This sectional unit is needed only if computeTransport 0 computeChemistry 0 or computeRND 0 integer dStep ut conc Steps between concentration outputs real T gt dt ut conc Time between concentration outputs integer nHistories t Number of 1 D time histories O 1 outFlags t i Monitor transport progress 66 O 1 outFlags_t_2 Monitor transport time step cuts O 1 outFlags t 3 Monitor transport dispersion solver integer outFlags t 4 Debug transport CMM trace back points O 1 outFlags t b Debug write transport conc O 1 outFlags t 6 Debug transport call and set up O 1 outFlags t 7 Debug transport CMM trace back integrals integer 0 4 outFlags t 8 Debug chemistry oll 1 outFlags t 9 Deb
51. N P T KEENAN M F WHEELER AND I YoTOV The application of mixed methods to subsurface simulation in Modeling and Computation in En vironmental Sciences R Helmig et al eds vol 59 of Notes on Numerical Fluid Mechanics Braunschweig 1997 Vieweg Publ pp 1 13 Enhanced cell centered finite differences for elliptic equations on general geometry SIAM J Sci Comput 18 1997 To appear T ARBOGAST C N Dawson D MOORE F SAAF C S SOUCIE M F WHEELER AND I Yotov Validation of the pics transport code tech rep Department of Computational and Applied Mathematics Rice University Houston Texas 1993 T ARBOGAST C N DAWSON AND M F WHEELER A parallel multiphase numerical model for subsurface contaminant transport with biodegradation kinetics in Computational Methods in Water Resources X Vol 2 A Peters et al eds Dordrecht The Netherlands 1994 Kluwer Academic Publishers pp 1499 1506 A parallel algorithm for two phase multicomponent contaminant transport Applications of Math 40 1995 pp 163 174 T ARBOGAST P T KEENAN M F WHEELER AND I YOTOV Logically rectangular mixed methods for darcy flow on general geometry in Proceedings of the 13th SPE Symposium on Reservoir Simulation held in San Antonio Texas February 12 15 1995 pp 51 59 SPE 29099 T ARBOGAST AND M F WHEELER A parallel numerical model for subsurface contaminant transport with biodegradation kinetics in
52. Similar formulas can be derived if several multi species phases are present As stated above the conversion assumes dilute solutions If the actual molar phase density changes drastically the computed equilibrium will differ from that given by 2 1 Certainly for the aqueous phase the approximation is almost always an excellent one Note that the solvent water must always be present in this formulation and it is checked that the water component has been specified Finally note that if desired the user can specify chmPhaseDensity 1 0 for a multi species phase to perform the calculation in mole fractions consider e g 2 5 Other units conversions will need to be made carefully elsewhere in the input data file s see 4 7 1 for the units conversion capabilities of the code 2 5 1 3 Conversion to equivalent reference chemical potentials When the conditions of chemical equilibrium are formulated as a constrained minimization problem it is the reference chemical potentials u IRS rather than the equilibrium constants that are the needed in the algorithms The conversion of the Nr equilibrium constants K1 Knp into a set of Ng equivalent reference chemical potentials is described in 28 The standard conditions of equilibrium result in the relation AG oi 2 8 RT 28 where the Standard Gibbs free energy change of the ith reaction has the form K exp Nc AGo V 10 ux ti Do Ag 2 9 j l The set
53. User s Guide to Parssiml The Parallel Subsurface Simulator Single Phase Version Parssiml v 2 1 May 1998 Minor updates to February 21 2006 Todd Arbogast Code Management The Center for Subsurface Modeling Mary F Wheeler director Steven Bryant associate director Todd Arbogast and Clint Dawson Texas Institute for Computational and Applied Mathematics The University of Texas at Austin Austin Texas 78712 mfw ticam utexas edu sbryant ticam utexas edu arbogast ticam utexas edu clint ticam utexas edu Copyright No portion of this document or program may be reproduced transmitted or otherwise copied without prior written permission of the Center for Subsurface Modeling CSM Texas Institute for Computational and Applied Mathematics The University of Texas at Austin except for the internal use of the CSM The authors make no representations or warranties about the correctness of any portion of the program code supplementary codes or associated documentation nor about the suitability of this program for any purpose and is in no way liable for any damages resulting from its use or misuse Table of Contents 1 Brief Summary 5 2 Governing Equations and Numerical Methods 6 2L Notation a 3o GP eae tok A Me ee ea ee ee el CUT oe eee ee 6 Did LO WES Ond he ot et Re eth sun SOND a hte AN eei o e c M 6 2 2 1 Viscosity Quarter Power Mixing Rule aaa aaa 7 2 22 Numerical Solution 234 ge igexue49 e Pee Be UR eh ROS Ge
54. ap constant 3 bcRegion z min Boundary region map 3 constant bcRegion_z_max Boundary region map SS Region 1 SSS Flow 2 bcType f BC type potential bcPresType in Dirichlet pressure potential head or hydraulicHead 1 nBCInterp f Number of interpolation times 1 0 bc f time BC time 4 atm bc f value BC value SSS Transport 1 bcType t BC type nBCInterp_t Number of interpolation times bc t time bc t value BC time BC conc value 1 0 0 2 2 ilyr O 1 yr 0 4 1 0 0 SS Region 2 SSS Flow 2 bcType f BC type 71 potential bcPresType_in Dirichlet pressure potential head or hydraulicHead 1 nBCInterp_f Number of interpolation times 0 bc f time BC time 3 atm bc f value BC value h SSS Transport O bcType t BC type SS Region 3 SSS Flow O bcType f BC type SSS Transport O bcType t BC type S WELLS 2 nWells Number of wells 1 maxnWelllInterp f Maximum number of flow interpolation times 2 t maxnWellInterp t Maximum number of transport interpolation times SS Well 1 1 wellType Well type 10 cm wellRadius Well radius 2 3 hWellIndex The x and y indices of well 3 zWellIndexLo Lower z index of well 4 zWellIndexHi Higher z index of well SSS Flow 1 nWellInterp f Number of interpolation times 1 0 wellInj_f_time Well time 10 USgal hr wellInj_f_value Well injection or extraction r
55. arameter The Runge Kutta Fehlberg convergence tolerance parameter used to control the estimated error lt integer gt rxn nDParams Number of real parameters for user s rxn fcns The number of real parameters needed to describe the specialized reactions 1 For each of the rxn nDParams real parameters give the following real 1 gt rxn dParams Real parameter The real double reaction parameters These are passed to the reaction functions in rxn c or rxnf f lt integer gt rxnnIParams Number of integer parameters for user s rxn fcns The number of integer parameters needed to describe the specialized reactions For each of the rxn_nIParams integer parameters give the following lt integer gt rxn_iParams Integer parameter The integer reaction parameters These are passed to the reaction functions in rxn c or rxnf f 5 12 Initial Conditions 8 INITIAL CONDITIONS 5 12 1 Flow SS Flow This sectional unit is needed only if computeFlow 0 lt real grid array L T gt tt velX X velocity A grid array of x Darcy velocity values given over the 3 D domain as point centered data in x and cell centered data in y and z for a total of nx 1 x ny nz velocities see 84 7 3 real grid array L T gt tt velY Y velocity A grid array of y Darcy velocity values given over the 3 D domain as point centered data in y and cell centered data in x and z for a total of nx x ny 1 nz velocities see 84 7 3
56. ate Y 888 Transport 1 nWelllInterp t Number of interpolation times welllnj t time welllInj t value Well time Well conc value 2 0 0 30 days 0 2 cc 1 0 0 1 0 0 SS Well 2 2 wellType Well type 15 cm wellRadius Well radius 3 4 hWellIndex The x and y indices of well 1 zWellIndexLo Lower z index of well 2 zWellIndexHi Higher z index of well SSS Flow 72 1 nWellInterp f Number of interpolation times 1 0 welllnj f time Well time 10 USgal hr wellInj f value Well injection or extraction rate h S LEAKS O nLeaks Number of leaks O maxnLeakInterp Maximum number of interpolation times S OUTPUT 1 runNumber Run number Sample data input file runDescription Title or description 1 verbosity Driver verbosity flag 0 debug Driver debug flag A outDir Output directory 1 outFormat3D Output 3 D file format 2 nSpeciesPerOutfile Number of species per 3 D output file 1 initial0ut Output 3 D initial conditions flag 1 outFormatiD Output 1 D file format SS Flow 1 dStep ut pres Steps between pressure outputs 1 dStep ut vel Steps between velocity outputs 1 dtOut_pres Time between pressure outputs 1 dt ut vel Time between velocity outputs potential presType_out potential pressure hydraulicHead or head 0 outFlags_f_1 Level of flow verbosity 0 outFlags_f_2 Level of flow debugging SS Transport 1 dStepOut_conc Steps betwe
57. ation R E Ewing ed no 1 in Frontiers in Applied Mathematics Society for Industrial and Applied Mathematics Philadelphia 1983 pp 35 106 Chapter II F SAAF A study of reactive transport phenomena in porous media PhD thesis Rice Univer sity Houston Texas 1996 5 A user s manual for nipsf Nonlinear interior point solver Fortran Tech Rep Tech Rep TR96 24 Department of Computational and Applied Mathematics Rice University Houston Texas 1996 F SAAF AND S BRYANT A user s manual for the geochemistry module in parsim1 reactive flow and transport code tech rep Center for Subsurface Modeling Texas Institute for Com putational and Applied Mathematics The University of Texas at Austin Austin Texas Jan 1997 W R SMITH AND R W MISSEN Chemical Reaction Equilibrium Analysis Wiley 1982 M F WHEELER T ARBOGAST S BRYANT C N DAWSON F SAAF AND C WANG New computational approaches for chemically reactive transport in porous media in Next Generation Environmental Models and Computational Methods NGEMCOM G Delic and M Wheeler eds Philadelphia 1997 Proceedings of the U S Environmental Protection Agency Workshop NGEMCOM SIAM pp 217 226 W I ZANGWILL Nonlinear Programming A Unified Approach Prentice Hall 1969 83 Index 21 21 21 21 21 21 19 20 E 18 18 20 21 52 c 18 37 e 18 18 20 21 52 h 18 18 i 18 18 25 t 18
58. cell of leak The leak occurs over a rectangular region Specify the three indices of the lower cell of the leak in the grid given as cell centered data indices integer integer integer leakIndexHi Upper cell of leak The leak occurs over a rectangular region Specify the three indices of the upper cell of the leak in the grid given as cell centered data indices For each of the nSpecies species give the following lt integer gt nLeakInterp Number of interpolation times The number of interpolation times used in the definition of the leak as a piecewise continuous linear function of time Constant extrapolation is used before and after the first and last times thus a single interpolation time gives a constant condition The times must not decrease however two pairs with the same time represent a discontinuous jump For each of the nLeakInterp interpolation times give the following lt real T gt leakInj t time Leak time The time at which the next interpolation value is to be used real 1 T gt leakInj_t_value Leak conc value The next interpolation value the molar rate of fluid leaking into the domain 47 5 16 Output S OUTPUT lt integer gt runNumber Run number The number of the current run Currently not used lt line of characters describing the run gt runDescription Title or description A title or description of the current run Maximum of 120 characters including spaces
59. ction well 45 45 relative permeability 33 relTol t 27 relTolIF f 27 p 6 RND see radionuclide decay RNSF 11 11 12 16 17 Robin condition 44 runDescription 48 Runge Kutta Fehlberg 10 41 runNumber 48 RXN see specialized chemistry rxn c 10 40 41 rxn comp 41 rxn dParams 41 rxn iParams 41 rxn nComps 41 41 rxn_nDParams 41 41 rxn nIParams 41 41 rxn nMicroSteps 41 rxn tolerance 41 rxnf f 10 40 41 scalar 33 scalar 3 34 section 20 sectional commands 19 20 25 52 sectional units see sectional commands 89 semicolon 19 serial 5 SF 11 11 12 15 16 gi 6 single species 36 SOLUTION PARAMETERS 26 sorp 6 8 36 Sorption 35 sorpType 37 source or sink 6 space 19 19 21 22 sparging well 45 specialized chemistry 5 9 10 10 26 40 SPECIALIZED REACTIONS 40 specie 37 specieName 37 specieNumber 40 stoich 12 38 38 stoichiometric formulation see SF stoichK 15 38 39 sub subsection 20 subdomain 5 27 28 subsection 20 subtraction 21 switchFlag 14 17 29 symmetric 33 symmetric 33 34 t 6 tab 19 19 21 tauMin 17 29 tdtMax_f 30 tdtMax t 31 Tecplot 5 48 56 57 temperature 21 26 36 Texas Institute for Computational and Applied Mathematics 5 tFinal 30 TIME 30 time 6 time splitting technique 6 8 9 17 tInitial 30 transDisp 8 35 transDispAry 8 35 Transport 27 80 42 44 46 49 transp
60. depend on the face considered e g for an x face the full cell centered range is 1 ny 1 nz A grid array that is to be fully specified is given as a data block see above A data block arising to specify a grid array can require a tremendous amount of information It may be convenient to put this information in an auxiliary data file in which case it could be read easily as INCLUDE lt file name gt It is easy to incorporate the auxiliary data file even if it has the wrong physical units Moreover it is easy to ignore any preamble that the file may contain using either the IGNORE_WORDS or IGNORE_LINES commands as in kg m 3 hr IGNORE LINES 2 INCLUDE file name gt which ignores the first two lines of the file lt file name and applies units conversion to the values Examples follow constant 2 7e 2 m nearly constant 2 7e 2 1 3 3 2 2 5 5 2 7e 2 1 2 cell or point 3 2 5 has half the value nearly constant 4 6 5 2 2 items per cell or point 2 1 3 2 4 5 5 2 7 2 3 2 3 4 5 3 1 6 overlapping areas have the values 3 1 and 6 23 nearly constant 2 5 1 1 2 0 3 6 7 a boundary grid array 2 5 304 7 800 3 5e2 fully specified 24 Chapter 5 The Input Data File s All input is read from data files that contain information arranged in sections subsections and sub subsections This data is read item by item and order is important Call the initial inpu
61. discussion on the meaning of the options lt integer gt nBCInterp f Number of interpolation times The number of interpolation times used in the definition of the boundary condition as a piecewise continuous linear function of time Constant extrapolation is used before and after the first and last times thus a single interpolation time gives a constant condition The times must not decrease however two pairs with the same time represent a discontinuous jump 43 For each of the nBCInterp_f interpolation times give the following lt real T gt bc f time BC time The time at which the next interpolation value is to be used lt real L T Neumann or M L T 2 Dir pres pot or L Dir head hHead gt bc f value BC value The next interpolation value 5 13 1 2 Transport SSS Transport This sectional unit is needed only if computeTransport 0 oll 111 2 bcType t BC type Flag declaring the type of the boundary condition for transport The command line argument t gives the options 0 for no flow condition or outflow condition where it should be the case that u v gt 0 and then DVc v 0 0 normal dispersive flux 1 for an inflow condition i e a Robin condition where it should be the case that u v lt 0 and then uc DVc v cgu v specified normal advective flux 2 for a Dirichlet condition c cg specified concentration Above cp needs to be specified as a function of time Function o
62. e file rxn c or written in Fortran in rxnf f 10 2 4 2 General Chemistry Package Reactive transport in porous media can be simulated with a full complement of homogeneous and heterogeneous reactions complexation adsorption ion exchange and precipitation dissolution of both equilibrium and kinetic type Kinetics are assumed to be of the mass action distance from equilibrium type where the reaction rate is proportional to the product of powers of component concentrations The package also permits generalised Monod style rate expressions This means that many biochemical biodegradation and bioremediation reactions can be incorporated directly into the general chemistry computation The chemistry reactions should be set up as if they were in a batch system that is the reaction rates should be given on a moles per unit pore volume i e flowing phase volume basis Reactions are not computed in the wells Three options exist for the choice of equilibrium module namely SF stoichiometric formula tion UNSF unreduced non stoichiometric formulation and RNSF reduced non stoichiometric formulation These options are described in detail in 25 As a general guideline we recommend the RNSF version as it is fastest for most problems that we have encountered The SF version has the advantage of exact mass balance and generalizes easier to multi phase problems but typically runs slower due to a larger number of variables in particular
63. e of the given file The error may be in a subsequently included file Error allocating memory on processor lt number gt There was an exhaustion of memory on the given processor You can first try to reduce the amount of workspace memory asked for If this fails you may be able to reduce the size of the problem you are solving Finally you can try to use a machine with more memory Error between input file lt file name gt line number lt number gt and input file lt file name gt line number lt number gt An otherwise unspecified error occurred somewhere between the given lines of the given files The error may be in an intervening included file Error in data block There is an error in a numerical data block see 4 7 2 on page 22 and also grid array in 4 7 3 on page 22 Perhaps the data block has the wrong number of items Error in input file lt file name gt before line number lt number gt An otherwise unspecified error occurred somewhere before the given line of the given file The error may be in a previously included file This message may be preceded by a more specific error message Error in input file lt file name gt on line number lt number gt An otherwise unspecified error occurred on the given line of the given file This message may be preceded by a more specific error message The item read was not of the appropriate type to continue e g an integer was expected but a real number was found The actual error
64. e the following lt real 1 gt stoichK Rate law powers Powers on the component concentration in the rate law 2 16 2 5 2 2 in the form ch Com ponents which appear in an expression of the form xi z are treated separately see next set of inputs and should be given a power of 0 here j For each of the nComps components give the following lt real 1 gt halfSatConst Half saturation constants The so called half saturation constants in Monod kinetics 2 16 2 5 2 2 which are the constants K in terms of the form Kt A component that does not appear in the rate expression in this form should be given a value of 0 here 5 9 Miscible Displacement S MISCIBLE DISPLACEMENT This sectional unit is needed only if computeFlow O and computeTransport 0 computeChemistry 0 or computeRND 0 oll 1 modelMiscDisp Miscible displacement model Flag declaring how the fluid viscosity will be computed in the flow computation The command line argument t gives the options 0 for a the fixed viscosity fluidViscosity 1 for the quarter power mixing rule 5 9 1 Quarter Power Mixing Rule Use the following when modelMiscDisp selects the quarter power mixing rulemodelMiscDisp This rule is defined above in 82 2 1 on page 7 For each of the nSpecies species give the following 39 lt real 1 gt viscosityRatio Viscosity ratio The ratio of the viscosity of the species to the resident
65. eactions can be extremely stiff and they are solved by a diagonalization of the system to increase stability and robustness The half lives of the species must be unique 17 Chapter 3 Program Name and Command Line Arguments The executable program for Parssiml is called simply parssim The unix calling specification is parssim i lt file name gt e lt file name gt E c lt file name gt t u v h The optional command line arguments are i lt file name gt The initial general input file lt file name gt or infile if omitted e lt file name gt Echo input to file lt file name gt or echo if omitted This is useful to track down errors in the input files It is also useful to strip most comments out of the data files and merge the included files into a single file since the echo file is in the form of an input file E Echo the input to standard out This is useful to track down errors in the input files c lt file name gt Echo general chemistry input data to file file name gt or echoChem if omitted This is useful to track down errors in the input description of the general chemical reactions if used t Just display Parssiml input type information These type numbers are used to input data u Just display Parssiml input physical units information This allows the user to see the predefined units available their dimensions and the initial internal base units
66. ecular weight 6 molecularWeight 6 37 Monod type kinetic 11 14 38 39 MPI 5 n 6 Ho 6 multi species 36 multi species 12 13 36 multiplication 21 nBCInterp f 43 44 nBCInterp t 44 44 nBCRegions 42 42 43 nChains 40 40 nComps 11 37 37 38 39 nDom_x 27 nDom_y 27 nDom_z 27 27 45 ndtMax_f 30 30 ndtMax_t 30 31 nearly constant 22 22 negation 21 Neumann condition 43 newline 19 19 20 21 nHistories_t 49 50 nHistSpecies t 50 50 nLeakInterp 47 47 nLeaks 47 47 no flow condition 43 44 non ideal 36 non ideal 36 nonlinear sorption 9 normal vector 6 notation 6 25 nPhases 36 36 nProducts 11 37 37 nSorpTypes 35 35 37 nSpecies 11 37 37 39 41 42 44 46 47 nSpeciesPer utfile 48 nStepsMax f 30 nStepsMax t 30 nStepsMax t per f 30 ntCutMax t 28 ntReact 17 29 v 6 numerical dispersion 10 nViol 16 29 nWellInterp f 46 46 nWellInterp t 46 46 nWells 45 45 nx 22 31 32 33 35 41 43 57 ny 22 31 32 33 35 41 43 57 nz 22 31 32 33 35 41 43 50 57 odeAlgType 17 29 order boundary regions 42 43 components 37 phases 36 38 products 37 sorption capacities 35 37 species 37 40 41 outDir 48 outFlags f 1 49 outFlags f 2 49 outFlags t 1 49 outFlags t 2 50 outFlags t 3 50 outFlags_t_4 50 outFlags_t_5 50 outFlags t 6 50 outFlags t 7 50 outFlags_t_8 50 outFlags t 9 50 outflow condition 43 44 outF
67. eded before INCLUDE file names S Equivalent to SECTION SECTION Begins a section of the data file Follow by the section name SS Equivalent to SUBSECTION 8SS Equivalent to SUBSUBSECTION SUBSECTION Begins a subsection of the data file Follow by the subsection name SUBSUBSECTION Begins a sub subsection of the data file Follow by the sub subsection name 4 3 Comments Explicit comments appear either after the comment symbol or between the BEGIN COMMENT and END COMMENT commands A comment begun with runs to the newline symbol so it can be used to remove lines or ends of lines from a data file without actually losing the information Note that 4t is an ordinary character if it is embedded in a string It is acceptable to nest BEGIN COMMENT and END COMMENT commands Implicit comments can be made before the first section Note that an unprocessed sectional unit is an implicit comment Moreover sectional units given bogus names such as remark can be used also as implicit comments 4 4 Sectional Units Sections subsections and sub subsections are initiated by a SECTION command SUBSECTION command or SUBSUBSECTION command respectively and a specified string It is important to note that entire sectional units may be skipped on input if they are not needed according to the options specified by the data This allows major routines to be omitted from the calculation such as flow transport
68. en concentration outputs 1 dt ut conc Time between concentration outputs 3 nHistories t Number of 1 D time histories 1 outFlags t 1 Monitor transport progress 1 outFlags t 2 Monitor transport time step cuts 1 outFlags t 3 Monitor transport dispersion solver O outFlags t 4 Debug transport CMM trace back points O outFlags t 5 Debug write transport conc O outFlags t 6 Debug transport call and set up O outFlags t 7 Debug transport CMM trace back integrals O outFlags t 8 Debug chemistry O outFlags t 9 Debug transport dispersion linear system assembly SSS History 1 histType t History type 1 CP 2 CT 3 MP 4 MT histIndexLo_t_x _y _z Lower cell of region histIndexHi_t_x _y _z Upper cell of region nHistSpecies_t Number of species in the history 3 2 histSpecie t Species in the history 1 dStepOut_hist_t Steps between history outputs 73 4 dtOut_hist_t Time between history outputs SSS History 2 histType t History type 1 CP 2 CT 3 MP 4 MT 2 4 histIndexLo t x _y _z Lower cell of region 2 4 histIndexHi t x _y _z Upper cell of region nHistSpecies t Number of species in the history 1 6 6 1 1 1 histSpecie t Species in the history 1 dStep0ut_hist_t Steps between history outputs i dtOut_hist_t Time between history outputs SSS History 3 1 histType t History type 1 CP 2 CT 3 MP 4 MT 2 2 4 histIndexLo_t_x _y _z Lower
69. es it 2 35 Transport As y a we v A Dix pem Bes ee eh de a on oa ee Ge 7 241 Diffusion Dispersion Tensor y Sk ee Be P Hoe ee bi DS e Reo 8 2 3 2 Linear Sorption Model e 8 2 3 3 Numerical Solution s uel tue A a ee 9 2 4 Chemistry Overview s o s ac u soea ello rss 10 2 4 1 Specialized Chemistry e 10 2 4 2 General Chemistry Package 2e 11 24 3 Numerical Solution 3E 2 ga a e dex a ee Ede xiu x EUR 11 2 5 General Ghemistry ose 4 4 oe geb i OR ee oe qom ee BU emu det PIE dos 11 2 5 bhermodynamics 4 2e ce Keka ple ple e liom RUE PR RIPE E D iid 11 2 5 2 Kinetic reactions and their rate laws 2 22 14 2 5 3 Numerical parameters oaoa 15 2 6 Radionuclide D c y 425 2 eso ok E oe RR UE RO REI 17 3 Program Name and Command Line Arguments 18 4 General Input File Syntax and Physical Units 19 4l Whitespace mu mL ek B E d osx Eo qo ee ERES 19 4 2 Commands i uxo a XU Rec gb S Rob E De Rod Re Ro ab ce ED RR a 19 4 3 COMMENDS zu E de has A A AIME EE Rc ee ae ee 20 AA Sectional Units m ep tU SEA eS Rue Eee ee ee Se QR EG E Ge BY 20 AD Jhist Delimiters os See Gi eae ce eR RD Ru RR a Pob eee le Re 20 4 6 KeyWords 21i Suis he hoe ee Ge PAV ele ee ep Ae Ma ae pss 21 AM Data Items zone mem sate ae Ue He ee a eee Sec ae RR oe a 21 4 71 Physical Uniis e da ao ee et Be he el ee ee e 21 4 1 2 Data Blocks and the Repetition Symbol rn 22 4 7 3 Gr
70. f the chain into this specie there are none for the first one for the second etc 5 11 Specialized Reactions S SPECIALIZED REACTIONS This sectional unit is needed only if computeChemistry lt 0 If computeChemistry gt 0 this section will be skipped The specialized reaction module is intended to be modifiable by the user so this section is quite generic The user provides nonequilibrium and nonequilibrium reaction functions in the C code file rxn c or in the Fortran code file rxnf f A single C preprocessor directive in rxn c must be set to select whether the reaction functions are written in C or Fortran This sectional unit uses the ordering of the species that is defined elsewhere 40 lt integer gt rxn nComps Number of component species in reactions The number of species that take part in the reactions This number may be less than nSpecies For each of the rxn_nComps component species give the following lt integer gt rxn_comp Species number of the next reaction component The species involved in the reactions are selected here and possibly reordered Order according to how the reaction functions in rxn c or rxnf f expect them Species are numbered beginning with 1 and ordered as defined elsewhere lt integer gt rxn nMicroSteps Number of reaction micro steps The number of sub steps to take in solving the reactions per transport time step lt double 1 gt rxn_tolerance Convergence tolerance p
71. f Time Specification The following lines appear only if a function of time needs to be specified that is bcType_t 0 For each of the nSpecies species give the following lt integer gt nBCInterp t Number of interpolation times The number of interpolation times used in the definition of the boundary condition as a piecewise continuous linear function of time Constant extrapolation is used before and after the first and last times thus a single interpolation time gives a constant condition The times must not decrease however two pairs with the same time represent a discontinuous jump For each of the nBCInterp_t interpolation times give the following lt real T gt i bct time BC time The time at which the next interpolation value is to be used real 1 L 3 bc t value BC conc value The next interpolation value 44 5 14 Well Specification S WELLS integer nWells Number of wells The number of wells Note that if there are wells nDom_z must be set to 1 above lt integer gt maxnWellInterp f Maximum number of flow interpolation times The maximum number of interpolation times used to specify the flow well conditions see below integer maxnWelllInterp t Maximum number of transport interpolation times The maximum number of interpolation times used to specify the transport well conditions see below 5 14 1 Single Well Description SS Well Include one such sectional unit
72. f information are given for larger values oll 1 outFlags t 9 Debug transport dispersion linear system assembly Boolean to select debugging output for transport diffusion dispersion linear system assembly 5 16 2 1 Concentration Histories SSS History Include one such sectional unit per history as specified by nHistories_t lt integer gt histType t History type 1 CP 2 CT 3 MP 4 MT The type of output desired The options are 1 for the pore volume average of primary phase concentration CP 2 for total primary phase moles CT 3 for the pore volume average of total moles MP 4 for the total moles MT lt integer gt lt integer gt lt integer gt histIndexLo t x jy z Lower cell of region The history sum or average occurs over a rectangular region Specify the three indices of the lower cell of the region in the grid given as cell centered data indices lt integer gt lt integer gt lt integer gt histIndexHi t x y z Upper cell of region The history sum or average occurs over a rectangular region Specify the three indices of the upper cell of the region in the grid given as cell centered data indices integer nHistSpecies t Number of species in the history The number of species to include in the concentration history For each of the nHistSpecies_t species give the following lt integer gt histSpecie t Species in the history 50 lt integer gt dStep ut hist t Steps between hi
73. fluid viscosity fluidViscosity po Use a zero viscosity ratio to indicate a species that contributes to the mass fraction computa tions but not to the viscosity directly i e the ratio is infinity Ignore completely any species with a negative viscosity ratio 5 10 Radionuclide Decay S RADIONUCLIDE DECAY This sectional unit is needed only if computeRND 0 lt integer gt nChains Number of RND chains The number of distinct radionuclide decay chains SS Chain Include one such sectional unit per chain as specified by nChains This sectional unit uses the ordering of the species that is defined elsewhere lt integer gt chainLen Length of chain The number of species in the chain linear branching chainType linear or branching State whether the chain is a simple linear chain each species decays into only one other species or a branching chain each species can decompose into all succeeding species For each of the chainLen species give the following lt integer gt specieNumber Specie number Give the next species in the chain This makes use of the order of the species given elsewhere lt real T gt halfLife Half life The half life stable is indicated by a negative half life The half lives of the species in the chain must be unique lt omit gt real 1 gt 4 branchRatios Branch ratios If chainType is branching give the branching ratio of the first second third member o
74. for the next component oll 1 1 2 reactionType Reaction type Flag to select the class of reaction The command line argument t gives the options 0 for equilibrium controlled 1 for mass action type kinetic 2 for Monod type kinetic Equilibrium Reaction Use this syntax if reactionType specifies an equilibrium reaction real 1 gt pK Log 10 equilibrium constant The base 10 logarithm of the equilibrium constant Mass action Type Kinetic Reaction Use this syntax if reactionType specifies a mass action type kinetic reaction real 1 gt real 1 gt pKf pKb Log 10 forward and backward rate constants For Kinetic reactions give the base 10 logarithm of the forward and backward rate constants 1 For each of the nComps components give the following real 1 gt stoichK Rate law powers Powers on the component concentration in the rate law If taken equal to stoich the clas sical form of the rate law results 38 Monod Type Kinetic Reaction Use this syntax if reactionType specifies a monod type kinetic reaction lt real 1 gt lt real 1 99 gt pKf pKb Log 10 forward amp backward rate constants For Kinetic reactions give the base 10 logarithm of the forward and backward rate constants The backward rate constant is included for generality it is often not used in biochemical pathways and a small value e g 99 should be assigned For each of the nComps components giv
75. for z velocity data conc for concentration data of the th component or a range of components There may be an appropriate graphics package extension appended to the file name Only a single vector velocity file or 3 component velocity files will appear as is appropriate for the graphics package The number is sequential and corresponds to the respective step number and time entry in the file presTime velTime or concTime for pressure velocity or concentration data respectively The grid will be stored in the file grid 000 if the graphics package supports a separate grid file otherwise the grid is stored in each individual output data file See the documentation on Eye and Tecplot for information on the format of these files Both of these formats represent all data on the same grid that consists of the centers of the cells not the grid points the centers of the boundary faces the centers of the boundary edges and the eight corners of the domain Thus the output grid covers the same domain as the input grid however the input and output grids contain different points Pressure and concentration data is cell centered data so in the output file each cell center data value is associated to the appropriate point of the output grid in the interior of the domain On the boundary of the domain the output grid points that are face centers of the input grid have associated with them the data value of the adjacent cell The
76. he arrays give values for the x faces y faces and z faces and are of size nx 1 ny nz nx ny 1 nz and nx ny nz 1 respectively 5 5 2 Dispersion SS Dispersion This sectional unit is needed only if computeTransport 0 The diffusion dispersion tensor is defined above in 2 3 1 on page 8 lt integer gt uniformDispersion Flag for uniform dispersion Flag for uniform dispersion if true i e 1 or nonuniform if false i e 0 Use this syntax if uniformDispersion is true i e 1 lt real L 2 T gt molDiff Molecular diffusion The molecular diffusion coefficient lt real L gt longDisp Longitudinal dispersion The longitudinal dispersion coefficient lt real L gt transDisp Transverse dispersion The transverse dispersion coefficient Use this syntax if uniformDispersion is false i e 0 real grid array L 2 T molDiffAry Molecular diffusion array The molecular diffusion coefficient array real grid array L gt longDispAry Longitudinal dispersion array The longitudinal dispersion coefficient array real grid array L gt transDispAry Transverse dispersion array The transverse dispersion coefficient array 5 5 3 Linear Sorption Retardation SS Sorption This sectional unit is needed only if computeTransport 0 computeChemistry 0 or computeRND 0 The linear sorption model is defined above in 2 3 and 2 3 2 see pages 7 and 8 lt integer gt nSorpTypes
77. he scheme has a CFL time step constraint so relatively small time steps must be taken The scheme has very little numerical dispersion and is locally mass and volume conservative Each time step though small is computationally relatively inexpensive 2 3 3 4 First Order Godunov Method FOG The higher order Godunov Method can be used without the postprocessing step resulting in the first order Godunov method for advection 2 3 3 5 Diffusion Dispersion This subproblem is solved by the cell centered finite difference technique used for the flow problem as described above A simple Jacobi preconditioned conjugate gradient technique is used to solve the resulting linear system 2 4 Chemistry Overview Chemistry can be simulated with either user supplied specialized reaction functions or the general chemistry package 25 27 2 4 1 Specialized Chemistry Specialized chemistry reaction functions can be supplied by the user compiled into the code and invoked The input specification consists almost exclusively of two arrays of real numbers and integers that can be used as parameters in the reaction functions The nonequilibrium reactions are solved by a fourth fifth order Runge Kutta Fehlberg tech nique The equilibrium functions are assumed to be solved analytically or iteratively by the user Parssiml provides no support for nonlinear solution of these equations The user supplies the reaction functions R either written in C in th
78. hemical species or a single chemical species A multi species phase is typically the flowing phase single species phases are typically minerals real 1 L 3 omit chmPhaseDensity Molar density The molar density of the phase in moles per fluid phase volume This item occurs only if phaseType is multi species 5 7 General Chemistry Properties S CHEMISTRY This sectional unit is needed only if computeChemistry gt 0 lt real temperature reading gt absTemp Equilibrium temperature The temperature at which equilibrium is computed See 4 7 1 for temperature readings This value is currently not used real M LT 2 chemPres Equilibrium pressure A representative pressure at which equilibrium is computed This value is currently not used ideal non ideal ideal ideal or non ideal Indicate whether all phases of the solution are to be considered ideal or non ideal Currently only the ideal option is implemented 36 5 8 Chemical Species Properties S CHEMICAL SPECIES lt integer gt nComps Number of components The number of chemical species that are components of the general chemistry reactions lt integer gt nProducts Number of products The number of chemical species that are products of the general chemistry reactions nSpecies nComps nProducts For ease of exposition and internal to the code we let nSpecies be the sum of nComps and nProducts If computeChemistry lt 0 so general chemistry
79. id Arrays Constant Nearly Constant and Fully Specified 22 5 The Input Data File s 25 DL General Information co a ee pan p AR eA A Bos 25 5 2 Algorithm Solution Parameters 2 0 0 0 0 es 26 AS A atte Aube enero ates a aa eee Delano eae ae 27 95 2 2 Erans port idem Rose SAE Be ek AS be A A EG bud 27 Dion ERES 4x iue ied heheh we yk d SLM dp A Re ees 28 5 95 Lime COMO 0 sot IY Pot tad Ants ce ee hae a M aae sien 30 DIL UEIOWOG d eR oe Bao eno RR ee e AUR A ea WR OR go 30 9 39 27 Transport imus eee d gexem rt eee eee he ey Ge nS wR Geet 30 pide Spatial Grid ia uy Rege sce on dr AT AA VOR CIR NS Pow 31 5 4 40 Uniform Rectangular Grid aoaaa 31 5 4 2 Nonuniform Rectangular Grid a 32 5 43 An XY rectangular Z variable Grid enn 32 5 4 4 A Fully Logically Rectangular Grid o 32 5 5 Porous Medium Material Properties 2e 33 5 5 1 Permeabilities or Cand tetivities sued a EE De EERE RE REY 33 D202 DISPETSION a omm iade d nae ao eu sore G5 A is payee te Re 35 5 5 3 Linear Sorption Retardation uteri RR eh Be RI Ey det IRR 35 5 6 Phase Properties 2222 e loh rss 36 5 6 1 Phase Description ana a a e a a a g a E e eE a 36 5 7 General Chemistry Properties es 36 5 8 Chemical Species Properties ee 37 5 8 1 Species Description a 2424 8 0444 4 8 n ba Ro BOE E RU POR ee 37 5 9 Miscible Displacement 4 ll less 39 5 9 1 Q
80. in 82 1 on page 6 A template of the input data file is also available it consists of the appropriate lines those in typewriter font below 5 1 General Information 8 GENERAL INFO oll 1 computeFlow Compute flow Boolean for whether a flow field will be computed If no flow field is computed the flow field should be given in 5 12 on initial conditions 25 01111121153 computeTransport Compute transport Flag declaring how the transport advection will be computed The command line argument t gives the options 0 for no transport 1 for the Characteristics Mixed Method CMM 2 for the higher order Godunov Method HOG 3 for the first order Godunov Method FOG for diffusion and dispersion only i e diffusion but not transport is computed it is rec ommended for consistency that one specify computeFlow 0 above a zero velocity field in 5 12 on initial conditions and either noflow conditions in 5 13 on boundary conditions or periodic boundary conditions in 85 4 0111112113 Il 1 computeChemistry Compute chemistry Flag declaring whether chemistry will be computed The command line argument t gives the options The full set of options is 0 for no chemistry 1 for the general chemistry RNSF option 2 for the general chemistry SF option 3 for the general chemistry UNSF option 1 for the specialized chemistry routines If option 1 2 or 3 is chosen the code must have been compiled for that option o
81. ints 8 MATERIAL PROPERTIES 9 8 m sec 2 gravity Gravitational constant constant 1 0 porosity Porosity SS Permeability scalar IGNORE_LINES 2 INCLUDE perm dat SS Dispersion 1 uniformDispersion Flag for nonuniform dispersion 0 0 molDiff Molecular diffusion 75 0 0 longDisp Longitudinal dispersion 0 0 transDisp Transverse dispersion SS Sorption 1 porosity sorp Sorption capacities for each type S PHASE PROPERTIES 1 2037e 8 fluidViscosity Flowing phase viscosity 1 0 fluidDensity Flowing phase mass density 4 nPhases Number of phases SS Phase 1 Aqueous phaseNames Phase name multi species phaseType multi species or single species 55 55 chmPhaseDensity Molar density SS Phase 2 Minerali phaseNames Phase name single species phaseType multi species or single species SS Phase 3 Mineral2 phaseNames Phase name single species phaseType multi species or single species SS Phase 4 Mineral3 phaseNames Phase name single species phaseType multi species or single species S CHEMISTRY 298 0 absTemp Equilibrium temperature 1 0 chemPres Equilibrium pressure ideal ideal ideal or non ideal S CHEMICAL SPECIES 4 nComps Number of components 3 nProducts Number of products SS Component 1 A specieName Specie name 100 0 molecularWeight Molecular weight 0 0 phaseDist The alpha parameter in Henry s Law SSS Chemistry 1
82. ion Typical bioreactions are irreversible and the backward rate constant k is usually set to zero 2 5 2 3 Kinetic reaction products It may happen that several kinetic reactions yield a common product or by product e g carbon dioxide or that a single reaction produces multiple products The framework described above makes no explicit accommodation of these possibilities but it can in fact be extended to these situations The trick is to define some common products or by products as components and assign them negative stoichiometric coefficients stoichK in the reactions for the other product species 2 5 3 Numerical parameters We now discuss appropriate settings for the chemistry solution parameters in 5 2 3 For some of these parameters useful default values are given there The parameter chmLoadBalFlag specifies whether automatic parallel load balancing of the computational work among processors should be attempted This creates parallel communication overhead but generally speeds up the computation in parallel 2 5 3 1 Interpretation of the results Let us begin by considering termination criteria for the equilibrium calculation The parameter chmMaxIter specifies the largest number of iterations allowed to satisfy the stopping criterion if this is number is insufficient an error message will result and the calculation fails Regardless of what version is used the iterates are required to satisfy lt chmAbsTol 2 17
83. is needed only if computeTransport 0 and wellType does not select an extractor or a bottomhole producer For each of the nSpecies species give the following lt integer gt nWelllnterp t Number of interpolation times 1 65 For each of the nWelllnterp t interpolation times give the following real T wellInj_t_time Well time real 1 L 3 gt wellInj_t_value Well conc value S LEAKS lt integer gt nLeaks Number of leaks lt integer gt maxnLeakInterp Maximum number of interpolation times SS Leak Include one such sectional unit per leak as specified by nLeaks lt integer gt lt integer gt integer leakIndexLo Lower cell of leak lt integer gt lt integer gt lt integer gt leakIndexHi Upper cell of leak For each of the nSpecies species give the following lt integer gt nLeakInterp Number of interpolation times For each of the nLeakInterp interpolation times give the following lt real T gt leakInj t time Leak time real 1 T gt leakInj t value Leak conc value h S OUTPUT lt integer gt runNumber Run number lt line of characters describing the run gt runDescription Title or description O 1 verbosity Driver verbosity flag oll 111 211 3 I 4 amp debug Driver debug flag word of characters outDir Output directory O 1 4 outFormat3D Output 3 D file format integer nSpeciesPer utfile Number of species per 3
84. is not used nSpecies is the num ber that is important and it can be given by dividing arbitrarily the species into components and products above 5 8 1 Species Description SS Component SS H Component SS H20 Component SS Product Include one such sectional unit per component as specified by nComps Follow by one such sectional unit per product as specified by nProducts General chemistry requires knowledge of whether the species is a fundamental chemical com ponent or a product of a reaction and sometimes requires knowledge of where the hydrogen or water component is This sectional unit induces an ordering of the components products and species that is used elsewhere Note that the species order places all components before all products that is the component order index is the same as the species order index and the product order index plus nComps is the species order index Species are numbered beginning with 1 lt word of characters gt specieName Specie name The name of the specie Maximum of 15 characters lt real M gt molecularWeight Molecular weight The molecular weight of the specie in mass units per mole This is used only for computing radionuclide decay and miscible displacement lt real 1 gt phaseDist The alpha parameter in Henry s Law The Henry s Law constant for linear sorption The linear sorption model is defined above in 2 3 and 82 3 2 see pages 7 and 8 This is the parameter a
85. ispersion solver outFlags t 4 Debug transport CMM trace back points outFlags t 5 Debug write transport conc outFlags t 6 Debug transport call and set up outFlags t 7 Debug transport CMM trace back integrals outFlags t 8 Debug chemistry outFlags t 9 Debug transport dispersion linear system assembly Oo0o0O0OO0O0Oo0oOorRrro dk HH HH GE GR OH HF A 3 2 Auxilliary file perm dat A 10 layered medium 40 40 1 1600 1 1600 3 79 1600 5 1600 1 1600 5 1600 24 1600 3 1600 24 1600 5 1600 24 A 3 3 Auxilliary file AB_IC dat 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 0000000000000000000000000000000000000000 0011111111111111111111111111111111111100 0011111111111111111111111111111111111100 0011111111111111111111111111111111111100 0011111111111111111111111111111111111100 Q O 1 1 L 111 171 1 1 1111114 11 114121410 117 114734111011 100 0011111111111111111111111111111111111100 Q0 1 1 141 111 11 11 1141411145 1 11111141 1111141111900 0011111111111111111111111111111111111100 0011111111111111111111111111111111111100 0011111111111111111111111111111111111100 001111111111111111
86. l processor Each subdomain is given roughly the same number of cells Each processor is responsible for the simulation only in its subdomain The individual processors send information to each other during the computation The program consists of the four main parts driver flow transport and chemistry The driver routines are responsible for the user interface input and output as well as managing the coupling between the flow and transport routines Chemistry is called from within the transport routines The flow is simulated with a package called Parcel 14 It allows for the simulation of incom pressible single phase saturated flow with wells on geometrically general domains but logically rectangular and it uses a locally conservative cell centered finite difference scheme The transport routine ParTrans allows the simulation of multiphase transport with linear sorption radionuclide decay simple specialized chemical reactions and wells on general ge ometry Transport is simulated using a locally conservative method of characteristics called the Characteristics Mixed Finite Element Method CMM or a Godunov Method The general chemistry routine handles both equilibrium and kinetic reactions For equilibrium reactions it uses an interior point algorithm for the minimization of the Gibbs free energy and is therefore relatively robust even when mineral phases precipitate into existence or dissolve away Output files can be written in
87. lt integer gt lt integer gt lt integer gt periodicBC x y z Boolean for periodicity Set to 1 if the grid covers a periodic domain in the given coordinate direction i e if there are periodic boundary conditons o Il 1 zIsDepth Direction of the z coordinate Boolean for whether the z coordinate points down depth or up height that is whether z increases with depth or with height We remark that internal to the code z points down uniform rectangular xy rectangular gridType uniform rectangular xy rectangular or The type of grid that is to be specified The grid must be rectangular or logically rectangular The options are e uniform for uniformly rectangular i e each cell has the same dimensions e rectangular for nonuniformly rectangular i e each cell is rectangular but they may vary in size e xy rectangular for a grid that is nonuniform rectangular in x and y but fully variable in z e to specify a fully logically rectangular grid this is the beginning of a data block This must be followed by one of the following sets of input data items 5 4 1 Uniform Rectangular Grid Use the following when gridType is uniformgrid lt real L gt lt real L gt xMin xMax Minimal and maximal x points The minimal and maximal x points of the computational grid lt real L gt lt real L gt yMin yMax Minimal and maximal y points The minimal and maximal y points of the computational grid
88. mponent 37 H20 Component 37 head 43 49 49 History 50 hydraulicHead 43 49 49 ideal 36 INITIAL CONDITIONS 41 Leak 47 LEAKS 47 linear 40 MATERIAL PROPERTIES 33 MISCIBLE DISPLACEMENT 39 multi species 12 13 36 non ideal 36 OUTPUT 48 Permeability 33 Phase 36 PHASE PROPERTIES 36 porosity 36 potential 43 49 49 pressure 43 49 49 Product 37 RADIONUCLIDE DECAY 40 rectangular 31 32 Region 43 scalar 53 34 single species 36 SOLUTION PARAMETERS 26 Sorption 35 SPECIALIZED REACTIONS 40 symmetric 33 34 TIME 30 Transport 27 30 42 44 46 49 uniform 31 31 Well 45 WELLS 45 xx 34 xy 34 xy rectangular 31 32 xz 34 yy 94 yz 34 zz 34 kinetic 5 38 KKT conditions 15 16 28 Leak 47 leak 47 leakIndexHi 47 leakIndexLo 47 leakInj_t 47 leakInj t time 47 leakInj t value 47 LEAKS 47 line of characters 21 48 linear 40 linear 40 linear sorption 5 8 8 9 35 37 list delimiter 19 20 22 23 52 53 logically rectangular grid 5 longDisp 8 35 longDispAry 8 35 mass action type kinetic 11 38 MATERIAL PROPERTIES 33 87 maxIter_t 27 maxIterIF_f 27 maxnBCInterp f 42 maxnBCInterp t 42 maxnLeakInterp 47 maxnWellInterp f 45 maxnWellInterp t 45 MISCIBLE DISPLACEMENT 39 miscible displacement 7 37 39 81 mobile 8 mobile specie 7 mobility ratio 7 modelMiscDisp 39 39 molDiff 8 35 molDiffAry 8 35 mole fraction 11 12 mol
89. n x min Boundary region map A grid array of boundary region identification numbers from 1 to nBCRegions given over the 2 D boundary face for minimal x This is cell centered data one per cell face see 4 7 3 for a total of ny nz integers integer grid array bcRegion x max Boundary region map A grid array of boundary region identification numbers from 1 to nBCRegions given over the 2 D boundary face for maximal x This is cell centered data one per cell face see 84 7 3 for a total of ny nz integers The following lines appear only if periodicBC y is 0 integer grid array bcRegion y min Boundary region map A grid array of boundary region identification numbers from 1 to nBCRegions given over the 2 D boundary face for minimal y This is cell centered data one per cell face see 4 7 3 for a total of nx nz integers 42 lt integer grid array gt bcRegion y max Boundary region map A grid array of boundary region identification numbers from 1 to nBCRegions given over the 2 D boundary face for maximal y This is cell centered data one per cell face see 84 7 3 for a total of nx nz integers The following lines appear only if periodicBC z is 0 integer grid array bcRegionz min Boundary region map A grid array of boundary region identification numbers from 1 to nBCRegions given over the 2 D boundary face for minimal z This is cell centered data one per cell face see 4 7 3 for a total of nx
90. nded that at most one of these be used at a time For immobile components we have a much simpler equation P cq c1 n torso dea where the influence of the well is completely local Equilibrium reactions are also supported as the infinite limit of increasing forward and backward rate constants that maintain a fixed ratio and subject to local mass balance constraints This is actually solved by a time splitting technique so the transport terms and the reaction terms are treated independently of each other In the limit equilibrium reactions amount to algebraic con straints on the concentrations and they are computed before or after transporting and kinetically reacting all of the species 2 3 1 Diffusion Dispersion Tensor The diffusion dispersion tensor is D u dmoil ul diongE u duas 1 E u where I is the identity tensor E v is the tensor that projects onto the v direction and whose i j component is UjUj E v i om Ivi and if uniformDispersion is used then do is mo1Diff diong is longDisp and dirans is transDisp and otherwise dmo is molDiffAry diong is longDispAry and dirans is transDispAry 2 3 2 Linear Sorption Model A simple linear sorption model is available It consists of the and o terms in the transport equations above The o vary with space representing the variation of the rock properties from point to point The o is a relative strength factor for the given species
91. ng lt integer gt rxn comp Species number of the next reaction component lt integer gt rxn_nMicroSteps Number of reaction micro steps double 1 gt rxn tolerance Convergence tolerance parameter lt integer gt rxn_nDParams Number of real parameters for user s rxn fcns For each of the rxn_nDParams real parameters give the following lt real 1 gt rxn_dParams Real parameter lt integer gt rxn_nIParams Number of integer parameters for user s rxn fcns 63 For each of the rxn_nIParams integer parameters give the following lt integer gt rxn_iParams Integer parameter S INITIAL CONDITIONS SS Flow This sectional unit is needed only if computeFlow 0 real grid array L T gt velX X velocity real grid array L T gt velY Y velocity real grid array L T gt velZ Z velocity SS Transport This sectional unit is needed only if computeTransport 0 computeChemistry 0 or computeRND 0 For each of the nSpecies species give the following real grid array 1 L 3 conc Molar concentration S BOUNDARY CONDITIONS lt integer gt nBCRegions Number of boundary regions integer maxnBCInterp f Maximum number of flow interpolation times integer maxnBCInterp t Maximum number of transport interpolation times The following lines appear only if periodicBC x is 0 integer grid array bcRegion x min Boundary region map in
92. nit time It should be nonnegative 5 14 1 2 Transport SSS Transport This sectional unit is needed only if computeTransport 0 and wellType does not select an extractor or a bottomhole producer For each of the nSpecies species give the following lt integer gt nWellInterp t Number of interpolation times The number of interpolation times used in the definition of the boundary condition as a piecewise continuous linear function of time Constant extrapolation is used before and after the first and last times thus a single interpolation time gives a constant condition The times must not decrease however two pairs with the same time represent a discontinuous jump 1 For each of the nWellInterp t interpolation times give the following real T gt welllInj t time Well time The time at which the next interpolation value is to be used 46 real 1 L 3 wellInj t value Well conc value The next interpolation value the molar concentration of the injected fluid 5 15 Leaking Source Specification S LEAKS lt integer gt nLeaks Number of leaks The number of leaking sources lt integer gt maxnLeakInterp Maximum number of interpolation times The maximum number of interpolation times used to specify the leak conditions see below 5 15 1 Single Leak Description SS Leak Include one such sectional unit per leak as specified by nLeaks integer integer integer leakIndexLo Lower
93. of reference potentials corresponding to a given set of equilibrium constants is not unique and we may without loss of generality assign arbitrary values to No of the reference potentials The most convenient choice is clearly us 0 3905 NGS 2 10 which allows us to write 2 9 as AG were 2 11 By combining 2 8 and 2 11 the equivalent chemical reference potentials for the system can be expressed in the simple form UNe i RT log Ki i 1 Nr 2 12 13 2 5 2 Kinetic reactions and their rate laws We allow fairly general rate expressions for the N x kinetic reactions Expressions derived from classical mass action equilibrium are supported as are generalized Monod type kinetic expressions 2 5 2 1 Distance from equilibrium expressions The distance from equilibrium has long been used as a measure of the driving force for change in the present application for change in chemical composition Specifically we express the kinetic rates rE as the difference between a forward and a backward rate No rE kf ej keno i l NE 2 13 j l where kf and k are the forward and backward rate constants respectively and p IRNeXN is a matrix of powers on the components in the rate law An important special case is the situation p i e the powers on the component concentrations are equal to their stoichiometric coefficients in the product species This is the classical form of the rate law Note that
94. oints real L gt real L gt yMin yMax Minimal and maximal y points real L gt real L gt zMin zMax Minimal and maximal z points Use the following when gridType is rectangular For each of the nx 1 x grid points give the following real L xp The x grid points For each of the ny 1 y grid points give the following real L gt yp The y grid points For each of the nz 1 z grid points give the following real L gt zp The z grid points Use the following when gridType is xy rectangular For each of the nx 1 x grid points give the following real L gt xp The x grid points For each of the ny 1 y grid points give the following real L gt yp The y grid points For each of the nx 1 ny 1 nz 1 z grid points give the following real of a data block L gt zp The z grid points Use the following when gridType is Y For each of the nx 1 ny 1 nz 1 x grid points give real of a data block L gt xp The x grid points For each of the nx 1 ny 1 nz 1 y grid points give the following real of a data block L gt yp The y grid points For each of the nx 1 ny 1 nz 1 z grid points give the following real of a data block L gt zp The z grid points S MATERIAL PROPERTIES real L T 2 9 8 m sec 2 gt gravity Gravitational constant real grid array 1 gt porosity Porosity 60 SS Permeability SS Conduc
95. on 16 The variables tauMin and chmRho are directly related to the interior point algorithms Variable tauMin describes the smallest percentage of the movement to the non negativity boundary that is allowed It should be taken close to but strictly less than unity The parameter chmRho is the fraction by which the perturbation parameter in the interior point method see 26 is decreased as the solution is approached It is always less than unity Typically a small value gives the fastest convergence but too small a value may lead to the global convergence problem of sticking to the boundary see 19 Note that the settings of tauMin and chmRho are of importance only when inequality constraints are present in the formulation Thus if the RNSF version is used to solve a problem without mineral phases these settings are of no consequence 2 5 3 3 ODE integration control The parameters ntReact odeAlgType and switchFlag control the time integration of the ODEs governing kinetic reactions and are only used if such reactions are specified by reactionType The variable odeAlgType controls the type of explicit time stepping scheme used for the inte gration Currently the choices 0 1 2 are supported corresponding to the use of Forward Euler FE first order scheme second order Runge Kutta RK 2 and fourth order Runge Kutta RK 4 The parameter ntReact defines a target time step to be used in the kinetic integration It specifies
96. on real 1 gt real 1 99 gt pKf pKb Log 10 forward amp backward rate constants For each of the nComps components give the following lt real 1 gt stoichK Rate law powers For each of the nComps components give the following lt real 1 gt halfSatConst Half saturation constants S MISCIBLE DISPLACEMENT This sectional unit is needed only if computeFlow 0 and computeTransport 0 computeChemistry 0 or computeRND 0 O 1 modelMiscDisp Miscible displacement model Use the following when modelMiscDisp selects the quarter power mixing rule For each of the nSpecies species give the following real 1 gt viscosityRatio Viscosity ratio S RADIONUCLIDE DECAY This sectional unit is needed only if computeRND 0 lt integer gt nChains Number of RND chains SS Chain Include one such sectional unit per chain as specified by nChains lt integer gt chainLen Length of chain linear branching chainType linear or branching For each of the chainLen species give the following lt integer gt specieNumber Specie number real T halfLife Half life omit real 1 gt branchRatios Branch ratios S SPECIALIZED REACTIONS This sectional unit is needed only if computeChemistry lt 0 lt integer gt rxn_nComps Number of component species in reactions For each of the rxn_nComps component species give the followi
97. on parameter 0 chmScaleFlag Diagonal scaling 0 chmTestSolFlag Test second order sufficiency conditions dl chmGuessType initial guess 1 auto 0 transported 1 previous 0 chmInterpFlag Use interpolation in the back track line search 10 nViol Number of non monotonic line searches 1 0e 4 chmAlpha Alpha parameter 1 hessFlag Analytical Hessian 0 995 tauMin Movement to the boundary factor 1 0e 3 chmRho IP reduction factor of perturbation parameter 1 ntReact Number of reaction steps per transport step 2 odeAlgType ODE integration O RK1 1 RK2 2 RK4 0 switchFlag Species switching 10000 dWorkspace c Double workspace for chemistry 10000 iWorkspace_c Integer workspace for chemistry S TIME 0 tInitial Initial time 01 2 tFinal Final time SS Flow 100000 nStepsMax_f Maximum number of flow steps 110 0 0 1 tdtMax f dtMax f Time of onset Time step size SS Transport 100000 nStepsMax_t Maximum number of transport steps 1 nStepsMax_t_per_f Maximum number of transport steps per flow step 2 0 0 0 000624 1 0 0 005 tdtMax_t dtMax_t Time of onset Time step size S GRID 40 40 1 nx ny nz Number of grid cells 0 0 0 periodicBC x _y _z Boolean for periodicity 1 zIsDepth Direction of the z coordinate uniform 0 0 1 5 xMin xMax Minimal and maximal x points 0 0 1 0 yMin yMax Minimal and maximal y points 0 0 1 0 zMin zMax Minimal and maximal z po
98. onentOrder for each grid cell Oth erwise give three grid arrays of single item cell centered data one for each tensor component again as ordered by permComponentOrder 5 5 1 3 Symmetric Tensor Permeabilities Use the following when permType is symmetric cells components permGrouping Group by cells or by tensor components Give the entire permeability tensor for each grid cell or give successively a single component of the tensor for the entire grid XX yy Zz xy xz yz any permutation gt permComponent rder Order of tensor s components Declare the order of the permeability tensor component data lt 1 or 6 real grid arrays L 2 permeability or L T conductivity gt perm Symmetric tensor permeabilities One or six grid arrays of symmetric tensor permeability values given over the 3 D domain as cell centered data six or one number per cell see 84 7 3 If the permGrouping is cells give the six components of the tensor as ordered by permComponentOrder for each grid cell Otherwise give six grid arrays of single item cell centered data one for each tensor component again as ordered by permComponentOrder 34 5 5 1 4 Face Permeabilities Use the following when permType is face lt 3 real grid arrays L 2 permeability or L T conductivity gt perm Face permeabilities Three grid arrays of permeability values given over the 3 D domain as face centered data one per face see 4 7 3 T
99. ormat1D 48 57 outFormat3D 48 56 OUTPUT 48 output data concentration 49 56 88 concentration history 49 51 57 pressure 49 56 velocity 49 56 outward unit normal vector 6 overlap region 28 28 p 6 padx 28 28 pady 28 28 padz 28 28 parallel 5 Parcel 5 81 parentheses 21 21 parssim 18 Parssim1 5 12 18 81 ParTrans 5 81 pcTypelF_f 27 Peaceman correction 45 periodic boundary conditions 5 26 periodic boundary conditons 31 periodicBC x 31 42 periodicBC y 31 42 periodicBC z 31 43 perm 6 33 34 35 permComponentOrder 34 34 Permeability 33 permeability 6 33 permGrouping 34 34 permType 33 34 35 Phase 36 PHASE PROPERTIES 36 phaseDist 6 8 37 37 phaseldentity 12 38 phaseNames 36 phaseType 36 36 pK 38 pKb 38 39 pkf 38 39 point centered data 22 41 42 pore volume 6 11 33 porosity 33 36 36 potential 6 46 49 potential 43 49 49 precipitation 11 pres 56 pressure 6 36 pressure 6 43 49 49 presTime 56 presType_out 49 primary phase 50 57 Product 37 product 37 production well 45 program name 18 proper white space 19 22 Y 6 q 6 quarter power mixing rule 7 39 RADIONUCLIDE DECAY 40 radionuclide decay 5 6 8 9 17 26 37 40 rate law 14 17 38 39 raw data 5 48 56 57 reactions 8 reactionType 17 38 38 39 rectangular 31 32 reduced non stoichiometric form see RNSF Region 43 reinje
100. ort 5 5 6 20 81 transport equation 7 type information 18 typewriter font 25 typewriter font 11 25 u 6 uniform 31 31 uniformDispersion 8 35 units 13 18 21 21 22 23 25 26 42 52 54 unix 18 48 unreduced non stoichiometric form see UNSF UNSF 11 11 12 16 vel 56 velTime 56 velX 6 41 56 velx 56 velY 6 41 vely 56 velZ 6 42 velz 56 verbosity 48 version information 18 viscosity 6 viscosityRatio 40 visual rendering 56 volRefine_t_x 28 volRefine_t_y 28 volRefine t z 28 volTol 28 Wi 6 Well 45 well 5 11 27 45 46 bottomhole production 45 extraction 45 inactive 45 injection 45 production 45 reinjection 45 sparging 45 wellInj f 7 46 welliInj_f_time 46 wellInj f value 46 wellInj t 8 46 wellInj t time 46 welllInj t value 47 wellPres 46 wellRadius 45 WELLS 45 wellType 45 45 46 90 white space 19 19 23 word of characters 21 36 48 workspace 27 29 30 x 6 xMax 31 xMin 31 xp 32 xx 34 xy 34 xy rectangular 31 32 xz 34 y 6 yMax 31 yMin 31 yp 32 33 yy 34 yz 34 z 6 zIsDepth 31 45 zMax 31 zMin 31 zp 32 33 zWellIndexHi 45 zWellIndexLo 45 zz 34 91
101. ory of one or more species over a fixed region of space These files are formatted so that they can be viewed directly with the graphics package Tecplot as set by outFormat1D These 1 D file names begin with concHist have the number of the history and then two letters describing the type of data recorded More specifically the names are one of concHist CP for the pore volume average of primary phase concentration of species m 5 dn mnVn 5 Qn Va n n concHist CT for the pore volume average of total moles flowing and sorbed of species m S n T Onin Cnn al 5 Pn Va n concHist MP for total primary phase moles of species m 5 Qn Cmn Vn TL concHist MT for the total moles flowing and sorbed of species m S n n On Goo Vas n where the is the history number and Vn is the volume of the nth grid element of the fixed region A final i j k or i j k i j k gives the cell centered data index numbers of the fixed grid cell or range of cells histIndexLo_t_x histIndexLo t y histIndexLo_t_z and if different histIndexHi t x histIndexHi t y histIndexHi t z See the documentation on Tecplot for information on the format of these files 57 Appendix A Sample Input Files A 1 An input template GENERAL INPUT DATA FILE TEMPLATE FOR PARSSIM1 This is a template of the input data file for PARSSIM1 version 2 1 S GENERAL INFO O 1 computeFlow Compute flow O 111 112 3 computeTransport Com
102. other output grid points are associated to averages of nearby boundary data values This procedure is used so that no false maxima or minima are represented in the output files however this has the effect that any visual rendering may appear to be a bit too regular near the boundary The velocities are not strictly cell centered for example the x velocity velX is face centered in x 56 and cell centered in y and z Each cell centered value is obtained as the average of the two face centered values of the cell and the boundary values are either given or derived as in the case of cell centered data The Eye format is about as compact as possible The grid is stored only once in the grid file and each data file contains a single data item pressure velocity component or species concentration Thus the storage is minimized but the number of files is maximized The raw data format is as compact as possible No grid is produced and each data file contains a single data item pressure velocity component or species concentration The order follows the logically rectangular structure of the grid and is given with the x coordinate varying most rapidly then y and finally z as in Fortran All files contain cell centered data and are of size nx ny nz except the velocity files which contain face centered data e g the x velocity file is of size nx 1 ny nz 7 2 Data that is 1 D in time One can write 1 D files of the concentration hist
103. ow o 1 outFlags t 1 Monitor transport progress Boolean to select monitoring progress in the transport routines 49 oO Il 1 outFlags_t_2 Monitor transport time step cuts Boolean to select monitoring time step cuts in the transport routines oll 1 outFlags t 3 Monitor transport dispersion solver Boolean to select monitoring iterations of the transport diffusion dispersion linear solver lt integer gt outFlags_t_4 Debug transport CMM trace back points Flag to select debugging of the transport trace back points in the CMM advection scheme if selected by computeTransport The value 0 means no output a negative value gives a 3 D plot of points for the Eye visualization package and a positive integer k between 1 and nz gives a logically horizontal 2 D cross section at index level k for the Draw program oll 1 outFlags t 5 Debug write transport conc Boolean to select debugging output for transport conc values oll 1 outFlags t 6 Debug transport call and set up Boolean to select debugging output for transport call and set up information o 1 outFlags t 7 Debug transport CMM trace back integrals Boolean to select debugging output for the CMM transport trace back integrals if the CMM advection scheme is selected by computeTransport integer 0 4 gt outFlags t 8 Debug chemistry Debug transport reactions Flag to select debugging output for the chemistry routines Value 0 gives no information and increasing amounts o
104. pute transport oll 1 11 211 3 1 computeChemistry Compute chemistry O 1 4 computeRND Compute RND lt units expression L gt baseLength Internal length base units lt units expression M gt baseMass Internal mass base units lt units expression T gt baseTime Internal time base units degK degC degF baseTemperature degKl degC or degF S SOLUTION PARAMETERS lt integer 1 gt lt integer 1 gt lt integer 1 gt nDom_x y _Z Parallel subdomain divisions 1 auto select SS Flow This sectional unit is needed only if computeFlow 0 lt integer 100 maxIterIF f Interface maximum number of iterations real 1 1 0e 6 gt relTolIF f Interface relative tolerance real 1 absTolIF f Interface absolute tolerance ol 11 5 5 gt pcTypelF f Interface preconditioner integer 48000 dWorkspace f Double workspace for flow SS Transport This sectional unit is needed only if computeTransport 0 integer 100 maxIter t Dispersion maximum number of iterations real 1 1 0e 6 gt relTol t Dispersion relative tolerance integer 10000 dWorkspace t Double workspace for transport real 1 1 gt cflFactor Numerical CFL factor SSS CMM This sectional unit is needed only if computeTransport selects the CMM advection scheme integer 2 ntCutMax t Maximum number of time step cuts integer 2 integer 2 integer 2
105. r an error will result o Il 1 computeRND Compute RND Boolean to invoke the radionuclide decay routines lt units expression L gt baseLength Internal length base units Set the internal length base units This overrides the initial values that were compiled into the code The initial default values can be retained by using the value 1 lt units expression M gt baseMass Internal mass base units Set the internal mass base units This overrides the initial values that were compiled into the code these can be retained by using the value 1 lt units expression T gt baseTime Internal time base units Set the internal time base units This overrides the initial values that were compiled into the code these can be retained by using the value 1 degK degC degF baseTemperature degKl degC or degF Set the internal temperature base scale and units See 4 7 1 for temperature readings This overrides the initial values that were compiled into the code these can be retained by using the value 1 5 2 Algorithm Solution Parameters S SOLUTION PARAMETERS 26 lt integer 1 gt lt integer 1 gt lt integer 1 gt nDomx y Z Parallel subdomain divisions 1 auto select The number of processors or subdomains in each of the three coordinate directions Note that their product must equal the total number of processors used in the simulation An automatic selection in any direc
106. r flow Joe Eaton developed the miscible displacement option Transport ParTrans was developed primarily by Todd Arbogast Ashokkumar Chilakapati and Douglas Moore Other authors include Lawrence C Cowsar Clint N Dawson Frederic d Hennezel Mary F Wheeler and Ivan Yotov Chemistry The general geo chemistry package was developed by Fredrik Saaf 25 27 It has been modified by Steven Bryant and Joe Eaton The development of this code was supported in part by the U S Department of Energy DOE through the Partnership in Computational Science PICS Grand Challenge Program administered through the Center for Computational Sciences CCS at the Oak Ridge National Laboratory ORNL other grants from the DOE the U S National Science Foundation and the State of Texas Governor s Energy Office 8l Bibliography 1 12 13 T ARBOGAST S BRYANT C Dawson F SAAF C WANG AND M WHEELER Com putational methods for multiphase flow and reactive transport problems arising in subsurface contaminant remediation J Computational Appl Math 74 1996 pp 19 32 T ARBOGAST A CHILAKAPATI AND M F WHEELER A characteristic mixed method for contaminant transport and miscible displacement in Computational Methods in Water Resources IX Vol 1 Numerical Methods in Water Resources T F Russell et al eds Southampton U K 1992 Computational Mechanics Publications pp 77 84 T ARBOGAST C N DAWSO
107. r in Henry s Law SSS Chemistry 4 phaseIdentity Participating phase 0 0 1 0 0 0 1 0 stoich 0 Equilibrium 3 0103e 1 pK Log 10 equilibrium constant S MISCIBLE DISPLACEMENT 0 None S INITIAL CONDITIONS SS Transport INCLUDE AB_IC dat A constant 1 0 B constant 1 0e 6 C TT constant 1 0e 6 4 D INCLUDE AB_IC dat AB s INCLUDE AB_IC dat AC s INCLUDE AB_IC dat DB s S BOUNDARY CONDITIONS 3 nBCRegions Number of boundary regions 1 maxnBCInterp f Maximum number of flow interpolation times 3 maxnBCInterp t Maximum number of transport interpolation times constant 3 bcRegion_x_min Boundary region map constant 3 bcRegion_x_max Boundary region map constant 3 bcRegion_y_min Boundary region map constant 2 bcRegion_y_max Boundary region map constant 3 bcRegion_z_min Boundary region map constant 3 bcRegion_z_max Boundary region map SS Region 1 SSS Flow O Noflow SSS Transport 1 Inflow 110 0 0 25 A 110 0 0 0 Jj B 110 0 2 0 Ct 1 0 0 1 75 D 1 0 0 0 0 s AB s 1 0 0 0 0 AC s 1 0 0 0 0 s DB s SS Region 2 SSS Flow 1 Neumann 110 0 1 0 888 Transport 0 Noflow outflow SS Region 3 SSS Flow O Noflow 888 Transport 0 Noflow outflow S WELLS 1 nWells Number of wells 1 maxnWelllInterp f Maximum number of flow interpolation times 1 maxnWelllInterp t Maximum number of tran
108. r messages that may be encountered pertaining to input Specific error messages will not be described here but only the following relatively generic error messages Arithmetic operation error in evaluating expression A units expression has an arithmetic error of some kind Beginning of list symbol not read A user declared number of inputs is now required This list must begin with the appropriate symbol Beginning of SECTION command not found A new section was supposed to begin This indicates that the previous section had extra data that was not processed Beginning of SUBSECTION command not found A new subsection was supposed to begin This indicates that the previous subsection had extra data that was not processed Beginning of SUBSUBSECTION command not found A new sub subsection was supposed to begin This indicates that the previous sub subsection had extra data that was not processed Binary addition or subtraction in units expression A units expression cannot involve addition or subtraction Negation or positivization is allowed Override by enclosing the units expression in double square brackets as in 1 1 52 End of list symbol not found A user declared number of inputs was required This list must end with the appropriate symbol Perhaps there are too many items in the list Error after input file lt file name gt line number lt number gt An otherwise unspecified error occurred somewhere after the given lin
109. rt 100 maxIter_t Dispersion maximum number of iterations 1 0e 6 relTol_t Dispersion relative tolerance 10000 dWorkspace t Double workspace for transport 1 cflFactor Numerical CFL factor 67 SSS CMM 2 ntCutMax_t Maximum number of time step cuts 2 2 2 padx pady padz Number of cells for subdomain overlap pad 10 volTol Volume discrepancy tolerance 1 1 1 volRefine_t_x volRefine_t_y volRefine_t_z Trace back volume refinement SS Chemistry 0 chmLoadBalFlag Parallel chemistry load balancing 100 chmMaxIter Maximum number of nonlinear iterations 0001 chmAbsTol Absolute KKT tolerance 1 0e 16 chmEpsConc Minimal concentration parameter 0 chmScaleFlag Diagonal scaling 0 chmTestSolFlag Test second order sufficiency conditions zr chmGuessType initial guess 1 auto 0 transported i previous 0 chmInterpFlag Use interpolation in the back track line search 10 nViol Number of non monotonic line searches 1 1 hessFlag Analytical Hessian 0 8 tauMin Movement to the boundary factor 1 0e 2 chmRho IP reduction factor of perturbation parameter 2 ntReact Number of reaction steps per transport step 2 odeAlgType ODE integration 0 RK1 1 RK2 2 RK4 0 0e 4 chmAlpha Alpha parameter switchFlag Species switching 1000 dWorkspace_c Double workspace for chemistry 1000 iWorkspace_c Integer workspace for chemistry S TIME 0 tInitial Initial time 3
110. s each may have units The number of values needed is the number that the array requires per grid cell or point A grid array that consists of a few repeated sets of values i e an array that is nearly constant can be specified easily For a 3 D array of grid cell or point values the specification is nearly constant primary value s gt number not that or those values 22 lt low x index gt high x index low y index gt high y index gt lt low z index high z index lt value s gt j where the punctuation is not necessary since it is merely white space the curly brackets are list delimiters that count the number number not that or those values of exceptional patches and the number of values needed is the number that the array requires per grid cell or point The meaning is as follows Fill the array with the first set of values primary value s gt given Next for the range of indices specified by each set of six integers overwrite the array with the given set of values lt value s gt Subsequent exceptional patches may overwrite previous patches For a 2 D array of grid cell or point values on a boundary face the specification is nearly constant primary value s gt number not that or those values lt low first index gt lt high first index gt lt low second index gt lt high second index gt lt value s gt where the meaning is similar but now the ranges
111. s of the KKT conditions see 18 used in the stop ping criterion as discussed in 2 5 3 Caution must be exercised with species of small total concentration see 2 5 3 1 real moles L 3 1 0e 16 gt chmEpsConc Minimal concentration parameter A lower bound on concentrations in the minimization algorithm see 2 5 3 and 2 5 3 1 This munber must be much less than chmAbsTol 28 o 1 lt 0 chmScaleFlag Diagonal scaling Boolean to attempt improvement in the condition number of the Newton system by a diagonal scaling procedure o 1 lt 0 chmTestSolFlag Test second order sufficiency conditions Boolean to test second order sufficiency conditions of optimality at the computed solution see 82 5 3 1 1 O I 1 lt 1 gt chmGuessType initial guess 1 auto 0 transported 1 previous Flag to select the initial guess for the equilibration computation The command line argument t gives the options 0 for the transported concentration values 1 for the previous concentration values 1 for an automatic selection 0 for pure equilibrium 1 otherwise oO 1 lt 0 chmInterpFlag Use interpolation in the back track line search Boolean to select interpolation in the back tracking algorithm used in the line search step otherwise use simple reduction see 2 5 3 2 lt integer 10 gt nViol Number of non monotonic line searches A control for the back track line search algorithm If 0 the standard
112. se mass density fluidDensity ci is the sorption capacity sorp v is the pressure potential pressure defined as v p pgz 2 2 Flow The flow equation represents conservation of mass in incompressible single phase flow and it is V u q where the Darcy velocity is k u Vy H and where q is related to wellInj_f The viscosity p is either yy or it is given as a function of the species concentrations by some mixing rule This allows simulation of miscible displacement 2 2 1 Viscosity Quarter Power Mixing Rule This rule see 24 is NEN 4 ME En pol bu Mte i 0 where is the viscosity uy is the resident fluid viscosity of species 0 fluidViscosity r is the ratio of the viscosity of the i species to the resident fluid viscosity reciprocal of the mobility ratio and is the mass fraction of the ith species The mass fraction is Ei cwi p The resident mass fraction g is given by solving We remark that the model should be valid whether or not the resident fluid is one of the species transported in the simulation We note that the fluid is assumed incompressible so the density of the fluid is not affected by its composition i e we assume dilute solutions This can cause inconsistencies if heavy or light components are present in the fluid in high concentrations 2 2 2 Numerical Solution Flow is computed independently of transport with independent time step sizes Generally speaking for time le
113. section key not found lt list of words gt The subsection name from the given list was not found It is likely missing misspelled or out of order Subsubsection key not found lt list of words gt The sub subsection name from the given list was not found It is likely missing misspelled or out of order Too many nested include files maximum lt number gt allowed While there is no limit on the number of files that may be included there is a maximum number that may be nested Try to finish a file before calling for succeeding files or have the code recompiled with a greater number of nesting levels allowed Unit expression error A units expression has an otherwise unspecified error Unit expression has the wrong dimensions A units expression has the wrong physical dimensions Override by enclosing the units ex pression in double square brackets as in 2 Units expression improperly formed or unknown unit A units expression is improperly formed It is possible that an unknown unit was given and the expression became confused Units expression is not positive A units expression is not positive Override by enclosing the units expression in double square brackets as in 1 Units expression is too long A units expression is too long You must shorten it or work out the conversion factor yourself Unknown command The command name was misspelled or the desired command is not supported Unknown unit A units expres
114. sion has an unknown unit or perhaps it was misspelled The command line argument u can be used to obtain a list of supported units New unit names and values can be added to the code and then it can be recompiled otherwise the numerical value of a unit in the base units can be used together with the double square bracket notation Unmatched end comment command An END_COMMENT command appeared without a matching BEGIN_COMMENT command 54 Unsupported input option The requested option is not supported Try another option To determine the supported options use this manual or perhaps the command line argument t to display the type numbers of various options 55 Chapter 7 The Output Data Files The user may write output data files as described in the chapter on input Chapter 5 These files may be 3 D in space or 1 D in time The user may also write various other files for testing and debugging purposes as described in the chapter on input 7 1 Data that is 3 D in space One can write 3 D files of the pressure velocities and concentrations These files are formatted so that they can be viewed directly with the graphics package Eye or Tecplot or the raw data is given depending on the setting of outFormat3D These 3 D file names have the structure data number Here data is one of the following pres for pressure data vel for vector velocity data velx for r velocity data vely for y velocity data velz
115. sport interpolation times SS Well 1 1 injector 0 05 radius 3 3 x and y indices 1 low z 1 high z SSS Flow 1 0 0 1 5 SSS Transport 78 1 0 0 0 25 A 110 0 0 0 B 110 0 2 0 C 1 0 0 1 75 D 1 0 0 0 0 ABCs 1 0 0 0 0 AC s 1 0 0 0 0 s DB s S LEAKS O nLeaks Number of leaks O maxnLeakInterp Maximum number of interpolation times S OUTPUT 1 runNumber Run number Fredrik 2Dtest problem runDescription Title or description 1 verbosity Driver verbosity flag 1 debug Driver debug flag outDir Output directory 1 outFormat3D Output 3 D file format 10 nSpeciesPerOutfile Number of species per 3 D output file 1 initial0ut Output 3 D initial conditions flag 1 outFormatiD Output 1 D file format SS Flow 1 dStepOut_pres Steps between pressure outputs 1 dStep Out vel Steps between velocity outputs 1 dt ut pres Time between pressure outputs 1 dt ut vel Time between velocity outputs potential presType out potential pressure hydraulicHead or head O outFlags f 1 Level of flow verbosity O outFlags f 2 Level of flow debugging SS Transport 1 dStep ut conc Steps between concentration outputs 0045 dt ut conc Time between concentration outputs nHistories t Number of 1 D time histories outFlags t 1 Monitor transport progress outFlags t 2 Monitor transport time step cuts outFlags t 3 Monitor transport d
116. story outputs The number of transport steps between writing concentration history output data files The value 1 is interpreted as positive infinity i e generate no output due reaching a certain step number real T gt dtOut hist t Time between history outputs The transport time increment between writing concentration history output data files The value 1 is interpreted as positive infinity i e generate no output due reaching a certain time 51 Chapter 6 Input Error Messages Errors in the input data file s can occur for a number of reasons such as the sectional units are out of order key words or sectional unit names are mispelled the wrong type of data or incorrect physical dimensions are assigned an incorrect number of items appears in a delimited list data items are assigned out of order general syntax is incorrect It is easy to understand how the code parses the input file by using the command line argument e or E This will automatically follow the included files and it can be extremely helpful in tracking down errors in the input data file s When searching for errors keep in mind that the error could occur well before the error manifests itself to the code since it is possible for data if it is of the proper type to be assigned to the wrong variables for some time The actual assignment to variables can be determined by setting the debug flag to a high enough value There are a series of erro
117. t 2 Oil specieName Specie name 90 g molecularWeight Molecular weight 2 phaseDist The alpha parameter in Henry s Law 2 sorpType Sorption type SSS Chemistry 1 phaseIdentity Participating phase O compCharge Charge SS Product 3 Mineral specieName Specie name 94 g molecularWeight Molecular weight 1 phaseDist The alpha parameter in Henry s Law SSS Chemistry 2 phaseldentity Participating phase 4 1 stoich Stoichiometry formula number 1 reactionType Reaction type 3 30 5 0 pKf pKb Log 10 forward and backward rate constants 1 0 0 1 stoichK Rate law powers h S MISCIBLE DISPLACEMENT 1 modelMiscDisp Miscible displacement model 70 8 1 1 viscosityRatio Viscosity ratio h S INITIAL CONDITIONS SS Flow This sectional unit is needed only if computeFlow 0 constant 10 ft day velX X velocity constant 0 velY Y velocity constant 0 velZ Z velocity SS Transport constant 1 cc conc Molar concentration constant 2 cc constant 2 cc S BOUNDARY CONDITIONS 3 nBCRegions Number of boundary regions 1 maxnBCInterp f Maximum number of flow interpolation times 2 maxnBCInterp t Maximum number of transport interpolation times constant 1 bcRegion x min Boundary region map constant 2 bcRegion x max Boundary region map constant 3 bcRegion y min Boundary region map constant 3 bcRegion y max Boundary region m
118. t file infile or use the command line argument i to specify the name of the initial file Include commands can be used to separate data into additional files if this is convenient or necessary Below is an item by item description of the input file divided into its various sections Manual section headings and indented remarks are not part of the input file s To help distinguish the input file text from merely manual text the typewriter font is used for input text The previous section describes general syntactical features of the input reading routines that may be exploited quite generally such as declaring units and including comments these features will not be repeated below Data items are described literally or by type Literal optional forms of an item are separated by only one option should be selected Item types names are enclosed in angled brackets as in integer Real numbers have physical units the dimensions are given as in real M LT 2 where L refers to length M refers to mass T refers to time In a few cases a suggested value is given after a colon as in integer 100 Comments are optional of course A boolean is a variable that is the answer to a yes or no question value 1 or not 0 is yes value 0 is no A flag is variable that provides a selection from several options as specified by some integer Sample input data files appear in Appendix A All symbolic notation is defined above
119. teger grid array bcRegion x max Boundary region map The following lines appear only if periodicBC y is 0 integer grid array bcRegion y min Boundary region map integer grid array bcRegion y max Boundary region map The following lines appear only if periodicBC z is 0 integer grid array bcRegion z min Boundary region map integer grid array bcRegion z max Boundary region map SS Region Include one such sectional unit per boundary region as specified by nBCRegions SSS Flow This sectional unit is needed only if computeFlow 0 O 111 112 bcType f BC type The following lines appear only if a function of time needs to be specified that is bcType f 0 omit pressure potential head hydraulicHead bcPresType in Dirichlet pressure potential head or hydraulicHead integer nBCInterp f Number of interpolation times 1 For each of the nBCInterp f interpolation times give the following real T bc f time BC time real L T Neumann or M L T 2 Dir pres pot or L Dir head hHead gt bc f value BC value 64 SSS Transport This sectional unit is needed only if computeTransport 0 O 111 112 bcType t BC type The following lines appear only if a function of time needs to be specified that is bcType_t 0 For each of the nSpecies species give the following lt integer gt nBCInterp_t Number of interpolation times
120. the smallest number of reaction steps per transport step that we consider adequate As an example ntReact 1 means that the integration of the ODEs will proceed with a time step equal to that used for transport Internally to the chemistry routines this time step size will always be attempted at the beginning of a step however it may be cut during the step if non negativity violations are encountered We point out that the variable chmEpsConc plays the role of a user specified approximate zero concentration in the kinetic integration We do not allow any kinetic process to decrease a species concentration below epsConc If switchFlag 0 the sub division into equilibrium and kinetic reactions will remain the same throughout the simulation This is the normal mode of operation However with the setting switchFlag 1 the algorithms will attempt to switch kinetic reactions that appear to have gone to completion to the equilibrium subset The equivalent equilibrium constant defined in 2 14 is used for this purpose We point out that if species switching is used all rate laws should be in the classical form as noted above otherwise the results are unpredictable 2 6 Radionuclide Decay Radionuclide decay reactions can be handled by the chemistry routines or by a separate radionuclide decay routine that is supplied If used it amounts to a time splitting of the radionuclide decay reactions These exponential growth and decay r
121. tion can be obtained by specifying the number 1 Note that if there are wells nDom_z must be set to 1 5 2 1 Flow SS Flow This sectional unit is needed only if computeFlow 0 lt integer 100 gt maxIterIF f Interface maximum number of iterations Maximum number of iterations for solving the flow interface problem lt real 1 1 0e 6 gt relTolIF_f Interface relative tolerance Relative tolerance for solving the flow interface problem lt real 1 gt absTolIF f Interface absolute tolerance Absolute tolerance for solving the flow interface problem Its value depends on the given problem since it is unscaled A very small number or even 0 is appropriate in most cases oll 115 lt 5 gt pcTypeIF f Interface preconditioner Flag to select the preconditioner for solving the flow subdomain interface problem The command line argument t gives the options 0 for no preconditioner 1 for diagonal Jacobi preconditioning 5 for a balancing preconditioner lt integer 48000 gt dWorkspace_f Double workspace for flow The amount of Fortran double precision workspace to allocate for the flow routines The flow routines will abort if this value is too small In that case raise its value and rerun the code The amount of workspace needed depends on the subdomain grid size 5 2 2 Transport SS Transport This sectional unit is needed only if computeTransport 0 lt integer 100 gt maxIter_t Dispersion maximum
122. tion of the pressure variable The meaning is as follows potential for the flow potential v p pgz M LT 2 pressure for p M LT 2 hydraulicHead for H 4 pg L head for the head h p pg L Internal to the code the variable pressure is assumed to be potential lt integer 0 1000 gt outFlags f 1 Level of flow verbosity Flag to set the level of computation monitoring in the flow routines Level 5 gives a summary of the iterations used to solve the interface problem of the domain decomposition solver lt integer 0 1000 gt outFlags f 2 Level of flow debugging Flag to set the level of debugging output in the flow routines 5 16 2 Transport SS Transport This sectional unit is needed only if computeTransport 0 computeChemistry 0 or computeRND 0 lt integer gt dStep ut conc Steps between concentration outputs The number of transport steps between writing concentration output data files The value 1 is interpreted as positive infinity i e generate no output due reaching a certain step number real T gt dt ut conc Time between concentration outputs The transport time increment between writing concentration output data files The value 1 is interpreted as positive infinity i e generate no output due reaching a certain time integer nHistories_t Number of 1 D time histories The number of 1 D time output data concentration history files to write as described bel
123. tivity This sectional unit is needed only if computeFlow 0 scalar diagonal symmetric face permType scalar diagonal symmetric or face Use the following when permType is scalar real grid array L 2 permeability or L T conductivity gt perm Scalar permeabilities Use the following when permType is diagonal by cells by components permGrouping Group by grid cells or by tensor components xx yy zz any permutation gt permComponentOrder Order of tensor s components lt 1 or 3 real grid arrays L 2 permeability or L T conductivity gt perm Diagonal tensor permeabilities Use the following when permType is symmetric cells components permGrouping Group by cells or by tensor components XX yy Zz xy xz yz any permutation gt permComponentOrder Order of tensor s components lt 1 or 6 real grid arrays L 2 permeability or L T conductivity gt perm Symmetric tensor permeabilities Use the following when permType is face lt 3 real grid arrays L 2 permeability or L T conductivity gt perm Face permeabilities SS Dispersion This sectional unit is needed only if computeTransport 0 lt integer gt uniformDispersion Flag for uniform dispersion Use this syntax if uniformDispersion is true i e 1 real L 2 T molDiff Molecular diffusion lt real L gt longDisp Longitudinal dispersion lt real L gt transDisp
124. ts of the number lt real gt and causes the value to be converted into the set of base units used internally in the code The command line argument e or E displays the conversion computation for diagnostic purposes The base units can be set by the user see 5 1 below A list of supported units with their physical dimensions and initial base unit numerical values can be found by using the command line argument u The arithmetic operations supported are integral exponentiation or multiplication division negation and parentheses C and Exponentiation by an integer with no physical dimensions only is allowed We also allow addition and subtraction though these should not appear in a units expression The usual rules of precedence apply exponentiation is performed first from left to right then multiplication and division and finally addition and subtraction although parentheses can be used to alter this order A leading or is allowed a preceding 1 is assumed Physical units are optional the set of internal physical base units are assumed if none are given If units are given the physical dimensions must be proper for the given quantity or an error will result Dimension checking can be disabled by using double square bracket notation as in lt number gt lt units expression gt in which case addition and subtraction are allowed This also
125. uarter Power Mixing Rule 0 020000 ee ee eee 39 5 10 Radionuclide Decay de pos esses d eux ok Ros OR RISE RUE ee RE dees 40 5 11 Specialized Reactions i 162 9 ae ru Oh ie RU T qe gu WR CHR Nr EORR 40 5 12 Initial Conditions 2 54 dt bees tek BS Pies Sea lee he nee uei Al Debs MEOW PE 41 5 12 2 Wransport 4 5 9 koX Re a Ro d r9 OY ROS SUELE 42 5 13 Boundary Conditions 42 5 13 1 Boundary Region Specification 2 43 Hila Well Specification su uei SW EE Be Gd ke Ma wRORU SE SN RUE 45 5 14 1 Single Well Description lle 45 5 15 Leaking Source Specification o osoo a 47 5 15 1 Single Leak Description eh 4T 5162 Outputs atin tz tea e oe e deem cere Pe un SEO a Dite e Doe p cede E eT qud 48 HU MEDI IPRC 49 9 16 2 Transport Vi a fae eg Rs eS es RES VERSER 49 6 Input Error Messages 52 7 The Output Data Files 56 7 1 Data that is 3 D in space le s Sos s 56 7 2 Data that is 1 D intime 22e Ss s 57 A Sample Input Files 58 AL An input templates ue SRL EAS uda RORIS dde Surge oe OR 58 A 2 5anipleinput file d i e i ea Cem Be alg Se auo mele ie fes Pug d 67 A 3 Sampleinputfile2 4 usu kac pog oo RE RO 3 Ro d d RR des 74 A 3 1 The main input file sss ss 74 A32 Ausalliary files perm dat unto Bt deter Gu bok Gd EP bo ee BS 79 A 3 3 Auxilliary file AB IC dat Authors and Acknowledgments Bibliography Chapter 1
126. ug transport dispersion linear system assembly SSS History Include one such sectional unit per history as specified by nHistories t integer histType t History type 1 CP 2 CT 3 MP 4 MT integer integer integer histIndexLo t x Jy _z Lower cell of region integer integer integer histIndexHi t x _y _z Upper cell of region integer nHistSpecies t Number of species in the history For each of the nHistSpecies t species give the following integer histSpecie t Species in the history h integer dStep ut hist t Steps between history outputs real T dt ut hist t Time between history outputs A 2 Sample input file 1 SAMPLE GENERAL INPUT FILE FOR PARSSIM1 version 2 1 S GENERAL INFO 1 computeFlow Compute flow 2 computeTransport Compute transport 1 computeChemistry Compute chemistry O computeRND Compute RND cm baseLength Internal length base units g baseMass Internal mass base units nin baseTime Internal time base units degC baseTemperature degKl degC or degF S SOLUTION PARAMETERS 1 1 1 nDom x y _z Parallel subdomain divisions 1 auto select SS Flow 100 maxIterIF f Interface maximum number of iterations 1 0e 6 relTolIF f Interface relative tolerance 0 absTolIF_f Interface absolute tolerance 5 pcTypeIF_f Interface preconditioner 100000 dWorkspace f Double workspace for flow SS Transpo
127. vel m 1 we solve implicitly V ut gee ut E gym um i A logically rectangular cell centered finite difference procedure is used to discretize the equa tions 14 This method handles tensor k accurately 13 and non rectangular geometry by a map ping technique 4 see also 11 8 12 3 The Glowinski Wheeler 20 domain decomposition solution procedure is used to solve the re sulting equations in parallel This involves solving an interface problem The subdomain linear system is solved directly 2 3 Transport The transport equation represents conservation of moles of a given species 7 For mobile species Oc p E V cu DVc RY c 5 Q aioi RN cr 2 Cn dci where RC represents the chemical reactions RN the radionuclide decay terms a and o are the linear sorption terms phaseDist and sorp described below and q represents the wells as dei d 6 q Ci where q is the positive part of q i e an injection well q is the negative part of q i e an extraction well and 6 is the injected concentration of the well wellInj_t Leaking sources add a specified amount of moles to specified cells per unit time but they do not affect the flow field It should be noted that if the general chemistry package is used as described below the coefficient of RC ci Cn is actually aigi thus care must be exercised when using both the chemistry package and the linear sorption terms it is recomme

User`s Guide to Parssim1 - The Institute for Computational

Contents

Download Pdf Manuals

Related Search

Related Contents