Home

os3D/GIMRT

1. gt ig OS3D GIMRT Manual 17 Steefel amp Yabusaki level of preconditioning equal to 0 was used see the WATSOLV manual This could be remedied by either setting the level of preconditioning to 2 or 3 or by simply tightening up considerably the tolerance and using a level of 0 This second option appears to give the best results 3 8 Convergence Criteria Convergence is obtained where all of the function residuals have been reduced below a tolerance set by the parameters atol and rtol in the driver routine gimrt f and in the file newton f in the case of the code OS3D These parameters are located in data statements at the top of the named routines The tolerance on function residuals is computed from tolmaz atol rtol C 72 where C refers to the primary species concentrations If they are changed the code should be recom piled The functions residuals when multiplied by the timestep delt in order to compare to the tol erances atol and rtol have units of moles per m bulk volume The default setting for the absolute tolerance is 107 mole m7 which corresponds to 107 mole liter Convergence will also be obtained where either the absolute value of the log concentration correction or the linear concentration correction drops below 10 5 i e approaching the roundoff error of the computer 3 9 Updating Mineral Properties When convergence is achieved the mineral volumes and surface a
2. gt using data base master25 else open iunit5 file home steefel GIMRT database mastertemp data status old err 334 write gt using data base mastertemp endif One of the options above involves a non default data base or location specified by datal ne which is discussed in the section on input files below The source code comes with Makefiles Makegimrt and Makeos 3d which should be modified slightly depending on the particular operating system to be used Copy the file Makegimrt or Makeos3d over to Makefile and uncomment the lines appropriate for the particular system on which the software is being installed e g Sun DEC or IBM These options at this time involve only slight differences in the compilation flags and in the case of the IBM AIX operating system the command to invoke the FOR TRAN compiler The second file which may have to be modified before compilation is params inc in the case of the GIMRT code or params os inc in the case of the OS3D code This file contains the parameter state ments which dimension the various arrays in the codes The codes are set up so that at run time the di mensions of the arrays are checked to make sure they are sufficiently large If the user needs to increase the dimension of the arrays or wants to decrease them so that they are just large enough the following parameters can be changed e mcomp 1 plus the number of independent chemical components
3. 1 0 0 0 0 0 gibbsite 31 956 3 3 0 ht 1 0 al 3 3 0 h20 6 7965 78 00356 kaolinite 99 520 4 6 0 ht 2 0 al 3 2 0 sio2 adq 5 0 h2o0 5 1007 258 16044 chalcedony 22 688 1 1 0 sio2 aq 3 7281 60 08430 quartz 22 688 1 1 0 sio2 aq 3 9993 60 08430 sio2 am 29 0 1 1 0 sio2 aq 2 7136 60 08430 null O 1 0 0 0 0 0 0602 24465 0 3 1 0 h20 1 0 h 1 0 hceo3 7 8136 44 00980 null O 1 0 0 0 0 0 0 0 ES E aper Steefel amp Yabusaki 48 OS3D GIMRT Manual 13 ACKNOWLEDGMENTS The present versions of the codes GIMRT and OS3D were written for the Environmental Dynamics and Simulation task which is part of the Environmental Molecular Science Laboratory EMSL construction project at the Pacific Northwest Laboratory in Richland Washington The codes are at least partly based on acode originally written while the senior author was a doctoral student at Yale University studying with Professor Antonio Lasaga The present version of the code however owes much to our collaboration with Peter Lichtner while at the Universit t Bern in Switzerland and since then while he has been at the South west Research Institute He was responsible for writing the algorithm in the subroutine database f which reads the thermodynamic data base and rewrites the stoichiometric coefficients and log K s in t
4. e 15 De 10 0e 8 0e 8 0e 5 0e 7 e 6 WBRRRAWRARO Pp The input includes the following fields Steefel amp Yabusaki 1 e 10 1 0e 10 tee 6 e 9 qa S 1 e 6 goethite 1 4 kaolinite 1 4 quartz 1 5e 1 albite 1 5e 1 gt 3 16 4 602 10 0 goethite 1l e 05 0 001 1 e 5 mi rig 1 e 8 7 1 e 9 1 e 8 at 1 e 8 1 5 wt 3 16e 4 602 10 0 goethite e Card 15 NCHEM number of geochemically distinct initial and or boundary conditions to be included In what follows the subscript K refers to the geochemical condition number while I refers to the primary species e NAMC ITYPE K GUESS 1 K CTOT I K NCON LK NAMC e mee w te name of primary species which must match that specified in Card 8 ITYPE LK constraint to be used in initialization 1 mass balance using total concentration ctot 2 charge balance using the named species 3 mineral specified by ncon equilibrium constraint 4 pressure of gas bars fixed by ctot j i POE SO CO A lt Steefel amp Yabusaki 36 OS3D GIMRT Manual 5 species is a surface complex developed on the mineral specified by ncon In this special case guess and ctot taken on different meanings Guess in the case of a surface complex gives the site density on the mineral in units of moles surface hydroxyls
5. i e we visualize the internal domain as cubic Negative and positive geochemical condition numbers play a similar role to the roles they play in Card 16 except that the two codes behave somewhat differ ently here If one is using GIMRT then specifying a negative geochemical condition number implies a Dirichlet condition This is unnecessary where the flow is entering the domain along this face since in this case the code will automatically use the upstream exterior node in computing the advective flux across the face But where the advective flux is zero specifying a negative geochemical number turns the boundary condition along this face from a no flux condition into a Dirichlet condition This is not possible in OS3D at the present time because of the way the transport algorithm is set up however As noted above if Dirichlet conditions are desired in the case of OS3D internal nodes immediately adjacent to the boundary should be used with negative numbers This method also has the advantage that one may specify individual nodes rather than forcing an entire face to have a single type of boundary condition 16 INTERNAL spatially distinct zones OVER AND ABOVE default k condition jx jy jz lower left jx jy jz upper right jx x direction jy y direction jz z direction 2 111 3 11 0 1 11 el i 17 GEOCHEMICAL CONDITION USED AT BOUNDARIES in following order jx 1 jx nx jy 1 1 jy ny 1 jz 1 1 jz nz 2
6. n n Hf gt GF ont cp pr E 56 where B r Max o Min 2 2r 4 57 and 01 cr 1 1 6x 1 01 62541 Similar calculations are performed for negative velocities and the method is generalized for multiple di mensions The technique is third order accurate in smooth regions oscillation free along concentration fronts and does not exhibit overcompression of fronts in the presence of physical dispersion Steefel amp Yabusaki 14 OS3D GIMRT Manual 3 5 2 Diffusion and Dispersion After the advection calculation OS3D proceeds to solve for diffusion and dispersion OS3D assumes that the principal direction of flow is aligned with the grid i e D is a diagonal tensor A centered nu merical scheme second order accurate in space and time is applied to the diffusion equation The dif fusion dispersion equations are solved with the WATSOLV sparse matrix solver of VanderK waak et al 1995 The user should understand that parameterizations of hydrodynamic dispersion attempt to account for ad vection resulting from variations in the velocity field that cannot be resolved by the numerical grid In this standard representation based on Ficks law dispersion is driven by the gradient of concentration and parameterized by a linear proportionality constant 1 6 the diffusion coefficient This is a gross simpli fication of a complex and poorly understood phenome
7. Calculate secondary species concentrations Calculate total concentrations Calculate mineral reaction rates Calculate Jacobian matrix Solve system of equations Update concentrations until converged Sequential Iteration No Update mineral volume fractions Figure 3 Program structure for OS3D which uses time splitting of the transport and reaction step or a sequential iteration between transport and reaction concentration boundary condition Atthe boundary nodes where the flow passes out of the domain the flux is also assumed to be purely advective so that it passes out of the system with the downstream exterior node having no effect on the system Note that when the flow across a boundary is non zero we always assume that the dispersive and diffusive fluxes are equal to zero E er q Zy a AA lt gt lt lt Steefel amp Yabusaki 22 OS3D GIMRT Manual The second type of boundary condition is the no flux condition which is applied when there is no flow across the boundary In this case both the advective flux and the dispersive diffusive fluxes are equal to zero The last type often called the first type or a Dirichlet boundary condition is normally applied when the ad vective flux is equal to zero but there is a dispersive or diffusive flux across the boundary if the advective flux is non zero then the code automatically assumes a zero dispersive diffusive flux This condition is many case
8. The chief advantage of the global implicit approach is that unlike the TVD method one is not restricted rigorously to the Courant condition which requires that mass not be transported more than a single grid Steefel amp Yabusaki 8 OS3D GIMRT Manual cell in any one timestep Although this is an important criterion to observe in order to maintain accuracy in transient problems it is less important once the aqueous concentrations approach a quasi stationary state In simulations we have carried out we find that the operator splitting error is minimal when the Courant number is less than about 0 2 so even if one does not use an explicit transport step minimizing operator splitting error requires a small timestep If the interest in the problem is to model water rock interaction over geological periods of time then transport errors associated with the transient propagation of concentration fronts are less significant In these cases the global implicit approach can offer a real advantage because of the ability to take larger time steps once the system relaxes to a quasi stationary state In contrast stability requirements dictate that the TVD algorithm which uses an explicit or forward Euler method use a timestep smaller than a Courant number of approximately 1 2 3 1 Global Implicit Approach In the one step method employed by GIMRT the transport and reaction terms are solved simultaneously One approach would be to solve directly for
9. cation of the executable image e g GIMRT The user should also make sure that the bin directory is included in their search path generally in the cshrc or login file The search path is set with the command set path 7 RUNNING GIMRT and OS3D OS3D and GIMRT both read the user provided input file threedin on startup in addition to a file containing the thermodynamic data base The file threedin includes various control parameters and physical parameters which are necessary in order to run the code The file also includes information on the initial and boundary conditions the minerals and species to be included and the rate data to be used in the simulation The input is organized as a series of cards which include comment lines followed by the fields where the user inputs data Comment lines begin all cards which are numbered Input data begins on the second line after the last comment line i e comments and input are separated by one space 7 1 Control Parameters Control and physical parameters for the simulation are specified in Cards 1 through 7 of the input file threedin A typical example is given below THIS IS THE INPUT data TO GIMRT AND OS3D 1 title up to 80 characters Short Course Test Problem 2 database 1 blank space for default 3 tstep max delt delt initial delt ttol ijacob bpeg 0 005 1 e 06 0 01 0 5 d 03 4 jpor constantpor porosity vdissmax vpptmax 0 0 1 0 01 0 001 5
10. jtemp tinit tgrad OS3D GIMRT Manual 29 Steefel amp Yabusaki 0 25 0 0 0 6 master iswitch igamma h 2 2 7 nstop prtint i i 1 nstop 1 3 0 20 0 30 0 40 0 50 0 e Card 1 LTITLE Title of simulation up to 80 characters enclosed in single quotes e Card 2 DATA1 String up to 100 characters with location of a user provided thermodynamic data base If no location is provided i e the code reads then the code defaults to either master25 dataif the input temperature given below is 25 C ortomastertemp data if another temperature is given e Card 3 TSTEP DELT TTOL IJACOB BPEG TSTEP maximum time step in years allowed DELT initial time step in years TTOL parameter controlling the tolerance on increasing the size of the time step Making ttol smaller decreases the rate at which the time step is increased for the purpose of controlling either the time truncation error or for numerical stability IJACOB flag controlling update of Jacobian O update Jacobian every Newton iteration 1 lt when delt 1 2 tstep update Jacobian every 3 Newton iterations and when delt tstep update Jacobian every 6 iterations BPEG maximum change in molality in any one Newton iteration Reducing bpeg in creases the under relaxation of the Newton step and may sometimes prevent non convergence although the total number of iterations required may actually increase e Card 4 JPOR CONSTANTPOR
11. 13 nzoney nvy i dyy i i 1 nzoney 0 16 1 00 14 nzonez nvz i dzz i i 1 nzonez 0 5 0 5 e Card 12 NZONEX NVX D DXX NZONEX number of zones in X direction with different grid spacing Steefel amp Yabusaki 34 OS3D GIMRT Manual NVX D number of grid cells using the grid spacing DXX grid spacing e Card 13 NZONEY NVY DYY D NZONEY number of zones in Y direction with different grid spacing NVY D number of grid cells using the grid spacing DYY D grid spacing in Y direction e Card 14 NZONEZ NVZ I DZZ 1 NZONEZ number of zones in z direction with different grid spacing NVZ I number of grid cells using the grid spacing DZZ D grid spacing An example in which uneven spacing is used in the X direction and even spacing in the Y direction might be 12 nzonex nvx i xx i i 1 nzonex 2 16 0 5 16 2 0 13 nzoney nvy i sayy i il 1 nzoney 1 16 1 00 14 nzonez nvz i dzz i i 1 nzonez 0 5 0 5 In this case we have one zone in the X direction with spacings of 0 5 meters using 16 grid cells and a second zone again using 16 grid cells with a spacing of 2 meters The grid spacing in the Y direction is a uniform 1 meter 7 6 Geochemical Conditions In Card 15 the code reads information on the distinct geochemical conditions to be specified in the problem Without indicating yet where these conditions will be located they may be initial conditions or boundary conditio
12. 1471 VanderKwaak J E Forsyth P A MacQuarrie K T B and Sudicky E A 1995 WATSOLV Sparse Ma trix Iterative Solver Package Versions 1 01 Unpublished report Waterloo Centre for Groundwater Research University of Waterloo Waterloo Ontario Canada 23 p Van Zeggeren F and Storey S H 1970 The Computation of Chemical Equilibria Cambridge Univer sity Press Cambridge 176 p Walter A L Frind E O Blowes D W Ptacek C J and Molson J W 1994 Modeling of multicom ponent reactive transport in groundwater 1 Model development and evaluation Water Resources Res 30 3137 3148 Wolery T J Jackson K J Bourcier W L Bruton C J Viani B E Knauss K G and Delany J M 1990 Current status of the EQ3 6 software package for geochemical modeling In Chemical Modeling in Aqueous Systems II ed D C Melchior and R L Bassett ACS Symp Ser No 416 104 116 Zysset A F Stauffer and T Dracos 1994 Modeling of chemically reactive groundwater transport Water Resources Res 30 2217 2228 Steefel amp Yabusaki 52 OS3D GIMRT Manual 15 DESCRIPTION OF THE INPUT FILE THREEDIN e Card 1 LTITLE Title of simulation up to 80 characters enclosed in single quotes e Card 2 DATAI String up to 100 characters with location of a user provided thermodynamic data base If no location is provided i e the code reads 7 then the code defaults to either master25 data if the input temperature given below
13. A A A be ees a 48 4 REFERENCES 2 4 Sa pelos pet be Rat aces oR ae Gee EE OR AG ENS 49 15 DESCRIPTION OF THE INPUT FILE THREED N 52 List of Tables 1 Sample data base read by GIMRT and OS3D for a single temperature of 25 C 47 List of Figures 1 Discretization scheme for a one dimensional problem using integrated finite differences 11 2 Program structure for GIMRT which uses a global implicit or one step method used to solve the multi component reaction transport equation 20 3 Program structure for OS3D which uses time splitting of the transport and reaction step or a sequential iteration between transport and reaction 21 EN me be oe carves E gt OS3D GIMRT Manual 1 Steefel amp Yabusaki 1 INTRODUCTION OS3D GIMRT is anumerical software package for simulating multicomponent reactive transport in porous media The package consists of two principal components 1 the code OS3D Operator Splitting 3 Dimensional Reactive Transport which simulates reactive transport by either splitting the reaction and transport steps in time i e the classic time or operator splitting approach or by iterating sequentially be tween reactions and transport e g Yeh and Tripathi 1989 Walter et al 1994 and 2 the code GIMRT Global Implicit Multicomponent Reactive Transport which treats up to two dimensional reactive trans port with a one step or global implicit approa
14. Ae 48 lt gt 7 wn A lt n lt os lt eae oe gt E A s us Steefel amp Yabusaki 12 OS3D GIMRT Manual 3 4 Transport Algorithm in Global Implicit Approach The addition of an advection term makes the set of partial differential equations much more difficult to solve stably and accurately e g Press et al 1986 Patel et al 1988 although the set of equations can still be cast in the general form of Equation 44 One way to represent the first derivative which appears in the advection term is to use a central difference method OU _ Ug Uy where the subscript P indicates that the derivative is evaluated at the nodal point P and Az is 1 2 0x u The truncation error for the central difference representation is said to be of the order Ax be cause in the Taylor series expansion upon which it is based only quadratic terms and higher are neglected Lapidus and Pinder 1982 In contrast a forward difference representation given by U _Ug Up dale en has a truncation error of O x that is all terms in the Taylor series higher than first order are truncated Despite the smaller truncation error of the central difference method the method may not be numerically stable for advection dominant problems because of the non physical oscillations which can occur in the vicinity of sharp concentration fronts Patankar 1980 Daus and Frind 1985 Frind and Germain 1986
15. Ne submatrices corre sponding to the primary species at any grid point That is if we are writing the conservation equation for 102 ag at a particular point in space the equation f3i0 Will be differentiated with respect to all Ne primary species e g 5402 ag Alt K etc As an example consider a one dimensional system described by 4 primary species which is discretized using 4 grid points this is for illustrative purposes only since hopefully one would not attempt to solve differential equations using only 4 grid points The matrix A can be represented compactly as a series of 4 x 4 submatrices A such that A11 A32 0 0 A21 Az A23 0 0 A32 A33 Aza 0 0 Aas As A 71 Those readers with some familiarity with numerical methods will recognize that this is the structure of any finite difference formulation of a partial differential equation e g the temperature equation except that in this case the entries are submatrices rather than single coefficients If we reduced the system to a single chemical component then the forms of the temperature and reaction transport matrices would be identical The submatrices need not be filled entirely with non zero elements but in general they will be dense in systems with either a large number of minerals or a great deal of complexation The diagonal submatrices A11 A22 etc contain contributions from both the heterogeneous reaction term Ren and from the total soluble concentration
16. Patel et al 1988 It appears therefore that the numerical stability of a particular finite difference formulation for the first derivative depends on the relative magnitudes of the advection and dispersion diffusion terms The rela tive importance of the two can be described by the dimensionless grid Peclet number v zx Pegrid D 51 where Az refers to the grid spacing at any particular point in space Only for Pegria lt 2 is the central dif ference formulation unconditionally stable Daus and Frind 1985 Frind and Germain 1986 Patel et al 1988 The problem of stability can be cured by using a forward difference formulation written in terms of the grid point itself P and the upstream or upwind node Press et al 1986 It should be pointed out however that the upstream weighted formulation in addition to having a larger truncation error than the central difference method achieves its greater stability essentially by adding numerical dispersion i e over and above the physical dispersion present in the problem In a one dimensional problem how ever numerical dispersion presents a problem only in tracking transient concentration fronts The effect does not appear in steady state one dimensional systems Patankar 1980 Numerical dispersion can be reduced by refining the grid in the vicinity of a sharp concentration front To take advantage of the greater accuracy of the central difference method at small values of P
17. Sons New York 568 610 Helgeson H C and Kirkham D H 1974 Theoretical predictions of the thermodynamic behavior of aqueous electrolytes at high pressures and temperatures I Summary of the thermodynamic electrostatic properties of the solvent Amer J Sci 274 1089 1198 Helgeson H C Garrels R M and MacKenzie F T 1969 Evaluation of irreversible reactions in geo chemical processes involving minerals and aqueous solutions II Applications Geochim Cos mochim Acta 33 455 482 Hindmarsh A C 1977 Solution of block tridiagonal systems of linear algebraic equations UCID 30150 Hooyman G J 1961 On thermodynamic coupling of chemical reactions Proc Nat Acad Sci 47 1169 1173 Johnson J W Oelkers E H and Helgeson H C 1992 SUPCRT92 A software package for calculating the standard molal thermodynamic properties of minerals gases aqueous species and reactions from 1 to 5000 bars and 0 to 1000 C Computers in Geosciences in press ere A hea we e qt lt Steefel amp Yabusaki 50 OS3D GIMRT Manual Kirkner D J and Reeves H 1988 Multicomponent mass transport with homogeneous and heteroge neous chemical reactions Effect of chemistry on the choice of numerical algorithm 1 Theory Water Resources Res 24 1719 1729 Lapidus L and Pinder G F 1982 Numerical Solution of Partial Differential Equations in Science and Engineering John Wiley and Sons New York 677 p Lasaga
18. and GIMRT 5 1 Units Although there are some advantages to allowing the input of data in any self consistent set of units many of the parameters used in GIMRT and OS3D are given commonly in certain units so we have decided to require these units All units for input parameters e g velocities reaction rates etc are given in the description of the input file threedin All distances are in meters The code assumes that the reaction rates for minerals are in units of mol m s these are converted in the code to mol m Diffusion coefficients are in units of m s All other references to time time for plot files timesteps etc are in years Velocities for example are in units of m yr Concentrations for aqueous species are in units of mol kg water i e molality The partial pressure of gases are specified in bars 5 2 Boundary Conditions At the present time there are three possible boundary conditions which can be used in GIMRT and OS3D The most straightforward way to describe the possible boundary conditions in the codes is to actually begin with the third type of boundary also sometimes called a Danckwerts boundary conditon This is a boundary condition which requires continuity of flux across the boundary Both codes assume that the flux at these boundaries is purely advective i e no diffusion or dispersion occurs upstream Where the flow is into the system the flux at the first interior nodes downstream
19. base This will occur for example if one were to include redox reactions without listing O2 among the gases in Card 11 One would also get an error message if one lists a species for which the necessary basis species are missing in the input file It is also fairly common to have problems in obtaining convergence during the initialization process This usually occurs when a mineral constraint is specified and the initial guess for the primary species concen trations is off by many orders of magnitude This is common in particular with redox species At this writ ing the initialization procedure in both OS3D and GIMRT rely on a damped Newton Raphson method which in future developments should be accompanied by a somewhat more robust method where the ini tial guess may or may not be good or is unnecessary altogether In certain other cases the initialization may fail to converge due to nonsensical input as for example where the specified total concentrations and mineral constraints would make charge balancing impossible If convergence problems during the initial ization are encountered it is usually advisable to replace the charge balance and mineral constraints with constraints on the total concentrations so that one can gradually improve the initial guesses for the basis species gt gt 2 Steefel amp Yabusaki 44 OS3D GIMRT Manual If non convergence occurs as the code steps through time one ca
20. cautionary note here is that one must include either among the primary or secondary species a species which is used in writing the reaction in the EQ3 EQ6 data base Note that this database is provided with the code in a reformatted form For example if one wants to include aluminum in the simulation then the species Alt must appear among either the primary or secondary species since that species is used as the primary or basis species in the EQ3 EQ6 data base In the same way one must include CO y among the list of gases if carbon is to be present in the system If this is not done the code will not be able to find the reaction when searching the data base One should not however include any secondary species which cannot be written completely in terms of the basis species specified in the input file The simplest way to grasp these features is to check the database to see which species are used to write the reactions involving a species of interest The species H20 presents a special case since in many cases it is present in such abundance that its ac tivity can be taken as unity If one wants the activity of H20 to be 1 then it should be left out of the list of primary and secondary species The code will still balance the reactions properly as if H20 were present If on the other hand one wants to carry out a mass balance on H20 then it should be included in the list of primary species and it will be treated like any other species except
21. is assumed to be just uC ez the same as the advective flux at the upstream exterior node In this special case where the dispersive and diffusive flux are neglected the third type of boundary condition is the same as a Dirichlet or fixed IO te mensan ss 25 gee eee o ES lt Steefel amp Yabusaki 20 OS3D GIMRT Manual Specify initial and boundar conditions Calculate concentrations of secondary species Calculate total concentrations ve Calculate mineral water reaction rate Update Jacobian Yes Calculate derivatives of total concentrations and reaction rates with respect to primary species Assemble residuals f x and Jacobian matrix over entire grid olve linear set of equations for concentration corrections NY Check for convergence Yes Update mineral volumes and areas Transient temperature and flow field Yes Update temperature and flo Figure 2 Program structure for GIMRT which uses a global implicit or one step method used to solve the multi component reaction transport equation wren et lt lt lt 2 renee A A wom gt OS3D GIMRT Manual 21 Steefel amp Yabusaki Specify initial and boundary conditions Begin time stepping Advect all species using TVD algorithm Diffuse and or disperse species Sequential Iteration Yes Add source term Reaction at individual space points
22. its surface area would equal its initial surface area Amo when Vm 1 if the porosity were constant 4 PROGRAM STRUCTURE 4 1 Global Implicit The basic structure of the global implicit approach GIMRT is summarized in Figure 2 After specifying the initial and boundary conditions the code begins stepping through time Within each time step a se ries of Newton iterations is required because of the nonlinearity of the equations Each Newton iteration consists of calculating the function residuals the f s given by Equation 55 and the partial derivatives of these functions the Jacobian elements with respect to the unknown concentrations Equation 59 and then assembling and solving the linear set of Ne x M algebraic equations in order to obtain the correc tions to the component concentrations If the basis switching option is turned on then the code checks that the basis species chosen are those that form a linearly independent set and that are present in the greatest concentration If more concentrated species are present than those currently used for the primary species then a basis switch is carried out at the appropriate grid cell i e the reactions are rewritten in terms of a new set of primary species Unlike the operator splitting approach in OS3D the residual functions and the Jacobian matrix are computed over the entire spatial domain before solving a global system of linearized equations Once convergence is achieve
23. node should be This however is not implemented at the present time OS3D and GIMRT differ slightly in their application of a Dirichlet boundary condition where the advec tive flux across the boundary is zero In the case of GIMRT the user may specify that an exterior node be a fixed concentration node and that it affect the interior node i e there will be a non zero diffusive and or dispersive flux through the boundary In OS3D however the code always assumes the disper sive diffusive flux across the boundary is zero due to the way the transport algorithm is formulated If one wants to use a Dirichlet condition in the case of OS3D one does so by fixing the concentrations at the first node or nodes interior to the boundary This will have the same effect as allowing an exterior node influence the interior domain except that one can think of the boundary as having been moved inward by one grid cell The procedure for doing this is described in the section below explaining the input file threedin To summarize the boundary condition options in OS3D and GIMRT a non zero flow across the boundary means that the code will use the upstream exterior boundary concentrations to give a purely advective flux at the boundary and the code will allow a fluid packet to pass out of the system with no influence from a downstream exterior nodes Exterior nodes across a boundary through which there is no flow will not influence the interior domain unless th
24. of mineral m No Nz Tm sgn lo0g Qm Km Amkm JI 2 28 1 1 In the input file containing the mineral fluid rate data the user has a choice of specifying the activities of those species which will affect the far from equilibrium dissolution rate In the case where no species are given the far from equilibrium dissolution rate becomes zeroth order since Qm Km is a small number under these circumstances The user however may provide a dependence on any of the either primary or secondary species concentration in solution In addition multiple parallel reactions may be specified This allows one to easily incorporate for example both a pH independent and a pH dependent i e proton promoted rate of dissolution for a mineral which are then summed to give the total dissolution rate One can also model ligand promoted dissolution with this scheme by specifying that the far from equilibrium rate be proportional to the concentration of reactive adsorbed ligand Several examples are given below The temperature dependence of the reaction rate constant can be expressed reasonably well via an Arrhe nius equation Lasaga 1984 Since many rate constants are reported at 25 C it is more convenient to write the rate constant at some temperature as F i R E 298 23 where E is the activation energy k s is the rate constant at 25 C R is the gas constant and T is tem perature in the Kelvin
25. of the full reaction step and Citi are the indi vidual species concentrations at the end of the second 1 2 transport step i e at the end of the full reaction transport time step L refers to the spatial differential operator defined by L Ve u DV 36 This scheme is reported to be second order accurate in time Zysset et al 1994 Note that in this for mulation all of the individual species are transported while the reaction step is still formulated in terms of the total concentrations which are functions of the individual primary species i e the solution of Ne unknowns is required If basis switching is not used one could transport the total concentrations in Equa tion 33 rather than than the individual species concentrations and follow this with a reaction calculation and then another transport of the total concentrations Where basis switching is used as it often is in the OS3D code however it is easier to transport all of the aqueous species so that the algorithm can construct any basis set it chooses from the individually transported species The form of the equation solved in the reaction step is the same as in Equation 32 except that in the operator splitting case we do not have the spatial operator i e the operator becomes U E 37 Yeh and Tripathi 1989 and Valocchi and Malmstead 1989 pointed out that in some cases the classic time splitting approach described above can lead to operator splitting err
26. possible to figure out the quirks of these codes and also the art of handling the highly nonlinear problems which arise in reactive transport If the code fails during the initialization process it is usually because of an error in the input file threedin It is a wise policy to keep a template of the input file handy which has not been corrupted so the user may remind themselves of the structure A common pitfall is to simply have the wrong format in the input file The file is set up so that many of the fields are expected even where they are not used This minimizes to some extent the changes in the actual format So for example as described above the code will try to read 3 character strings containing the velocity fields even where the user has specified that either there is no advective flux or where the flow velocity is constant These files need not exist in this case however since the code will not try to open them unless the parameter iveloci ty described above is set to 1 Another common problem is a misspelling of a species name The user may check the databases mas ter25 data and mastertemp data and any other they have generated to verify that their spelling is correct If the message appears that a species has not been found in the data base then either the species has been misspelled or a species has been listed which is not accompanied somewhere else in the list by the appro priate species used in writing the reaction in the EQ3 data
27. precipitation rate by specifying dependences on the dissolution rate far from equilibrium Josep Soler of Yale University recently added the capability for differing precipitation and dissolution rate laws to an earlier version of the code so this capability may be added to next release of the code These expressions can also handle nonlinear expressions by adjusting the parameters n and M in Equation 28 By providing dependences on the dissolution rates far from equilibrium where most ex perimental studies to date have been carried out one can investigate a range of reaction mechanisms in cluding pH dependent and pH independent rates a dependence on the chlorinity of the solution and even ligand promoted dissolution i e it is possible to provide a dependence on a surface complex The rate law may depend on either secondary or primary species One can add other rate expressions however by altering the source code This involves adding or substituting the appropriate function to the overall rate for a mineral computed in the subroutine reaction f In addition one must add the derivative of this function with respect to the various basis species to the Jacobian calculation in the subroutine rctder f 12 4 Providing and or Updating a Velocity Field As described above in the sections describing the input file threedin the user may specify that a ve locity file should be read in Note that since the velocities are needed at the faces of grid
28. small amounts of a mineral are present than it is possible that dissolution of the entire amount present within a grid volume may oc cur in a single time step This typically results in either inaccurate propagation of a mineral front or in oscillations in the mineral volume fraction profile A k A Steefel amp Yabusaki 18 OS3D GIMRT Manual Initial reactive surface areas are specified along with initial mineral volume fractions by the user at startup Reactive surface areas of minerals are updated differently at the present time depending on whether the minerals are primary or secondary For primary minerals the surface areas are given by id H dissolution precipitation where m o is the initial volume fraction of the mineral m and o is the initial porosity of the medium This formulation ensures that as the volume fraction m of a mineral 0 its surface area also 0 In addition for both dissolving and precipitating minerals the term go 3 requires that the surface area of a mineral in contact with fluid gt 0 when the porosity of the medium 0 For secondary minerals i e those initially with a volume fraction 0 the reactive surface areas as a function of time are given by 2 3 y ra bm dissolution 2 3 76 precipitation m Am 0 which implies that if a secondary were to begin to dissolve
29. that 1 2 of the reaction rate is included as a source term in advance of the reaction module and so iteration between transport and reaction is required to achieve convergence Figure 3 Mineral volume fractions and the porosity are updated after all of the species have been reacted at all the space points in a similar fashion to the global implicit 5 FEATURES OF THE CODE Before running the code the user of OS3D and GIMRT needs to decide what physical process es he or she wants to model and design the simulation accordingly This involves choosing the boundary condi tions for the problem e g no flux boundaries or fixed concentration boundaries and setting up the initial grid and or initial mineral zones The user must decide how many grid points are needed to solve a par ticular problem Typically one proceeds by carrying out the simulations on a relatively coarse grid e g 30 to 50 grid points and only refining the grid later when one has a general idea of the behavior of the system The user also needs to decide what level of chemical complexity they wish to consider This can involve choices in the number of independent components primary species the number of secondary species and the number of reacting minerals The discussion in the sections which follow is designed to make the basis for these decisions clearer These sections are then followed by a step by step discussion of the input file threedin which is needed to run OS3D
30. the case of internal zones specified in Card 16 a TEO dr E E EE EIA 7 aggre x OS3D GIMRT Manual 37 Steefel amp Yabusaki positive geochemical condition number indicates that the rectangular domain which is specified should be initialized to the appropriate geochemical condition 2 in the example below but that the concentrations in this zone will not remain fixed once the transient calculation begins In contrast if a negative geo chemical condition is specified 2 in the example below the code will fix the concentrations to those specified by the appropriate geochemical condition number i e 2 in this case for the duration of the tran sient calculation As discussed in the section on Initial and Boundary Conditions this option is useful particularly when using the code OS3d if one wants to specify a Dirichlet boundary fixed concentration boundary We do so by fixing the first interior node in from the boundary thus effectively sliding the location of the boundary in one grid cell This is not necessary to do where the flow of fluid is into the interior domain at this point since in this case both OS3D and GIMRT will grab the upstream exterior node with the appropriate geochemical condition specified in Card 17 In Card 17 6 different condition numbers must be provided even in OD cases since the code will expect to read these The geochemical conditions are those to be specified along an entire face of the domain
31. the total concentration of surface hydroxyls on this mineral is set to a small number 107 10 mol kg and its concentration only builds up once the secondary mineral precipitates How each of the above initialization constraints is applied is described in the section on the input file threedin below 5 9 Temperature Gradients Since GIMRT and OS3D are based on finite difference methods they have no problem handling tem peratures which are non constant either in time or in space although this option has been implemented only in GIMRT due to time constraints on the preparation of the software If either a non default ther modynamic data base or the master25 data data base is not used at startup then the code defaults to the mastertemp data file where equilibrium constants are given at temperatures of 0 C 25 C 60 C 100 C 150 C 200 C 250 C and 300 C at pressures corresponding to the water saturation curve The code then fits a polynomial to the thermodynamic data so that equilibrium constants can be generated at any temperature between 0 C and 300 C along the water saturation curve Simulations are possible at other pressures and temperatures but the user must provide his or her own data base in the same format as master25 dataormastertemp data In addition the algorithm for fluid density must be modi fied as well Presently the fluid density is calculated as a function of temperature from a polynomial fit to the data in Helge
32. velocity components This field should include 3 character strings even if the code is to be run with a constant velocity otherwise a read error will result In the files since the velocities are face centered the X component of the velocity should have dimensions of 0 nx 1 ny 1 nz the Y component 1 nx 0 ny 1 nz and the Z component 1 nx 1 ny 0 nz e Card 21 ALFL ALFT longitudinal and transverse dispersivities in meters respectively e IDIFFUS DCOEFF PORPOW IDIFFUS flag controlling choice of a diffusion coefficient O default 1x 107 m s and an activation energy of 5 kcal mole 1 user provided uses dcoeff DCOEFF user provided diffusion coefficient Do in units of m s only used if idi f fus D PORPOW the cementation exponent in Archie s Law Equation 5 OS3D GIMRT Manual 41 Steefel amp Yabusaki 8 PARAMETERS IN FORTRAN FILES Some of the parameters which the user may want to change are located within the FORTRAN source files themselves These are normally parameters which may not require change as commonly as those in threedin Note however that if any of these parameters within the source files are changed the code must be recompiled with the make command The code GIMRT uses the WATSOLYV software package VanderKwaak et al 1995 to solve the global system of equations The code OS3D uses the package to solve the dispersion diffusion equation The various values of these para
33. where a user should place a routine for calculating or inputting a temperature are marked in the routine gimrt f 5 10 Updating Porosity and the Transport Coefficients A user can choose at runtime whether to include a feedback between the changes in the volume fractions of the various minerals in the system and the diffusion coefficient In the present version of the code there is no provision for a feedback between porosity and permeability This is planned as a feature of the next release To do so requires that one solve the continuity equation the equation for the conservation of fluid mass in addition to the reaction transport equation Steefel and Lasaga 1990 Steefel and Lasaga 1994 Local changes in the porosity result normally in local changes in the permeability thus causing a change in the distribution of the hydraulic head or fluid pressure The continuity equation should be solved in a one dimensional system with constant head boundary conditions and not constant flux conditions as was done in Berner 1980 The fluid velocity becomes infinite as the porosity gt O when the constant flux conditions are used If the user decides to update porosity during the course of the simulation then the magnitude of the dif fusion coefficient in the porous medium D changes according to Archie s Law Dullien 1979 D 78 where m is again the cementation exponent in Archie s Law This formulation ens
34. 1 1 1 1 1 16 INTERNAL spatially distinct zones IN PLACE OF default k condition jx jy jz lower left jx jy jz upper right jx x direction jy y direction jz z direction 2 111 111 0 1 11 1 1 1 17 GEOCHEMICAL CONDITION USED AT BOUNDARIES in following order Steefel amp Yabusaki 38 OS3D GIMRT Manual jx 1 jx nx jy 1 Jy ny 3251 jz nz 2 1 2 1 1 1 In the first example geochemical condition 2 is used along the left hand face of the domain i e JX 1 while geochemical condition 1 is used at all other faces However only those faces which have flow entering the domain from outside or those which specify a negative geochemical condition number none in the first example will actually influence the interior domain In the second example we fix grid cell 1 to be geochemical condition 2 for the entire duration of the transient calculation thus effectively turning it into a Dirichlet boundary condition at grid cell 1 This would be redundant however if the flow were coming into the domain at this point in the X direction since Card 17 tells the code to use geochemical condition 2 in computing the advective flux across the left hand face JX 1 Note that we have specified at the face corresponding to JY 1 that geochemical condition 2 to be used and that it be treated as a Dirichlet condition Only GIMRT will actually do anything with this If this option is tried in OS3D the code will r
35. A C 1981 Rate laws in chemical reactions In Kinetics of Geochemical Processes ed A C Lasaga and R J Kirkpatrick Rev Mineral 8 135 169 Lasaga A C 1984 Chemical kinetics of water rock interactions J Geophys Res 89 4009 4025 Lichtner P C 1985 Continuum model for simultaneous chemical reactions and mass transport in hy drothermal systems Geochim Cosmochim Acta 49 779 800 Lichtner P C 1988 The quasi stationary state approximation to coupled mass transport and fluid rock interaction in a porous medium Geochim Cosmochim Acta 52 143 165 Lichtner P C 1992 Time space continuum description of fluid rock interaction in permeable media J Geophys Res 28 3135 3155 i Marsily G de 1986 Quantitative Hydrogeology Acad Press New York 440 p Patankar S V 1980 Numerical Heat Transfer and Fluid Flow Hemisphere Publ New York 197 p Patel M K Cross M and Markatos N C 1988 An assessment of flow oriented schemes for reducing false diffusion Int Journ Numerical Methods in Engineering 26 p 2279 2304 Press W H Flannery B P Teukolsky S A and Vetterling W T 1986 Numerical Recipes The Art of Scientific Computing Cambridge University Press Cambridge 818 p Reed M H 1982 Calculation of multicomponent chemical equilibria and reaction processes in systems involving minerals gases and an aqueous phase Geochim Cosmochim Acta 46 513 528 Steefel C I 1992 Coupled fluid flow and c
36. IGAMMA MASTER string up to 18 characters enclosed in single quotes specifying master variable used in controlling time step ISWITCH flag controlling basis switching 0 no basis switching solves for logarithms of concentrations 1 not allowed at this time 2 basis switching every Newton iteration Improved method over earlier versions which involves solving for the linear concentrations and using re laxation to prevent overly large Newton steps or negative concentrations Results in reduced number of cases where non convergence occurs espe cially in redox examples Does require some extra CPU time however IGAMMA flag controlling method for computing activity coefficients 0 unit activity coefficients except for H20 2 extended Debye Hiickel without calculation of the derivatives of the ac tivity coefficients with respect to the primary species e Card 7 NSTOP PRTINT NSTOP number of graphics files to be printed out 0 no time stepping distribution of species only 21 number of output files at prtint intervals PRTINT times at which output should be written at least nstop times should be listed 7 2 Specifying Primary and Secondary Species In the next part of threedin Cards 8 through 9 the initial choice of primary and secondary species is specified Each list is terminated with x An example input would be 8 PRIMARY SPECIES co 2 tracer h OS3D GIMRT M
37. New York 764 p Bear J 1979 Hydraulics of Groundwater McGraw Hill New York 569 p Berner R A 1980 Early Diagenesis A Theoretical Approach Princeton Univ Press Princeton 241p Bowen R M 1968 On the stoichiometry of chemically reacting materials Arch Rational Mech Anal 29 114 124 Bruges E A Latto B and Ray A K 1966 New correlations and tables of the coefficient of viscosity of water and steam up to 1000 bar and 1000 C Internat Jour Heat and Mass Transfer 9 465 480 Daus A D and Frind E O 1985 An alternating direction Galerkin technique for simulation of contam inant transport in complex groundwater systems Water Resources Res 21 653 664 Dullien F A L 1979 Porous Media Academic Press San Diego 396 p Frind E O and Germain D 1986 Simulation of contaminant plumes with large dispersive contrast Evaluation of alternating direction Galerkin models Water Resources Res 22 1857 1873 Gupta A D Lake L W Pope G A Sephernoori K and King M J 1991 High resolution mono tonic schemes for reservoir fluid flow simulation Jn Situ 15 289 317 Helgeson H C 1968 Evaluation of irreversible reactions in geochemical processes involving minerals and aqueous solutions I Thermodynamic relations Geochim Cosmochim Acta 32 853 877 Helgeson H C 1979 Mass transfer among minerals and hydrothermal solutions In Geochemistry of Hydrothermal Ore Deposits ed H L Barnes John Wiley and
38. PNWVL 11 EL OS3D GIMRT Software for Modeling Multicomponent Multidimensional Reactive Transport USER MANUAL ao amp O M PROGRAMMER S GUIDE 4 x lt Sm C I Steefel S B Yabusaki Pacific Northwest National Laboratory Richland Washington 99352 U S A 1 Department of Geology University of South Florida Tampa Florida 33620 U S A E mail steefel margaux cas usf edu Copyright 1995 1996 by Battelle Memorial Institute All Rights Reserved Version 1 0 May 6 1996 DISCLAIMER This report was prepared as an account of work sponsored by an agency of the United States Government Neither the United States Government nor any agency thereof nor any of their employees make any warranty express or implied or assumes any legal liability or responsibility for the accuracy completeness or usefulness of any information apparatus product or process disclosed or represents that its use would not infringe privately owned rights Reference herein to any specific commercial product process or service by trade name trademark manufacturer or otherwise does not necessarily constitute or imply its endorsement recommendation or favoring by the United States Government or any agency thereof The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof fois te et gt wees Soha ea Ay DISCLAIMER P
39. VDISSMAX VPPTMAX JPOR flag controlling update of porosity and flow 0 no update of porosity 1 porosity is updated CONSTANTPOR porosity of the medium for the case where the porosity is not updated por 0 If jpor 0 the porosity is taken from constantpor rather than from the sum of the volume fractions of the minerals in the system If jpor 1 porosity updated constantpor is ignored VDISSMAX maximum rate of decrease of a mineral s volume fraction units of Can be used to limit the amount of mineral dissolution which can occur within a single time step Useful where either the dissolving phase is very soluble or where very small amounts are present Where very small amounts of a mineral are present inaccurate mineral front propagation may occur where most or all of the mineral present within a grid volume dissolves in a single time step Peeps mene e ra o wee gt se gt Steefel amp Yabusaki 30 OS3D GIMRT Manual VPPTMAX maximum rate of increase of a mineral s volume fraction units of yr e Card 5 JTEMP TINIT TGRAD JTEMP flag controlling temperature option O isothermal 1 linear temperature gradient operating only in X direction TINIT temperature in C at grid point 1 TGRAD temperature gradient in C m A positive gradient causes the temperature to increase from grid point 1 a negative gradient causes it to decrease e Card 6 MASTER ISWITCH
40. anual fe 3 al 3 8102 aq nat tal heo3 gt feoh tke 9 SECONDARY SPECIES oh aloh 2 al oh 2 al oh 3 aq al oh 4 0602 aq 1c03 2 h3sio4 nacl aq gt feo gt feoh2 gt feoh co kt e Card 8 PRIMARY SPECIES NAMC J name of primary species to be included enclosed in single quotes til marks end of primary species read 241 Steefel amp Yabusaki In this example H is chosen as the primary species implying that OH has been chosen as a secondary species A secondary species should be included only if it can be written completely in terms of the in cluded primary or basis species The exception to this is a species which is written in terms of H2O aq The reaction will be correctly balanced whether or not one chooses to explicitly include H20 among the primary species The inclusion of H20 ag only determines whether the code carries out a mass bal ance on H20 ag Since the reactions involving carbonate in the EQ3 data base are written locally in terms of CO g it should be included in the list of gases towards the end of the input file Similarly the inclu sion of redox reactions requires that Oz y be included in the list of gases since the EQ3 database writes many of the redox reactions using that species Species beginning with a gt as in gt feoh are surface complexes which need to b
41. ative geochemical condition number is used only by GIMRT transforming what would otherwise be a no flux boundary into a Dirichlet boundary see section on Initial and Boundary Conditions e Card 18 IPATH parameter only used when nxyz 1 i e no spatial discretization and controlling whether the mineral volume fractions are updated When ipath 1 andnxyz 1 the code runs in reaction path mode in order to model the movement of a fluid packet through the rock leaving reaction prod ucts behind When ipath 0 the code runs in true batch mode and the reaction products remain in the closed system e Card 19 PLOTDUMMY I character string for species to be plotted at a specific point in space e INTERVAL JXPLOT JYPLOT JZPLOT interval at which time dependent data the species listed above in plotdummy will be written to the output file breakthrough out followed by the location in X Y and Z space of the point of interest For a reaction path calculation jxplot jyplot and jzplot should be all equal to 1 r marks end of read e Card 20 IVELOCITY flag indicating whether to use a constant velocity field or to use the information found in the files specified below O assume constant velocities given in qxinit qyinit and qzinit below 1 use velocity components found in the files listed below AE Steefel amp Yabusaki 56 OS3D GIMRT Manual e QXINIT QYINIT QZINIT X Y and Z components of the cons
42. cells the X components of the velocity must go from 0 nx 1 ny 1 nz the Y component from 1 nx 0 ny 1 nz and the Z component from 1 nx 1 ny 0 nz The velocities are to be stored so that they can be read with an implicit do loop of the form for a 3D problem read 23 qx jx jy jz jx 0 nx jy 1 ny jz 1 nz read 24 qy jx Jy J z jx 1 nx jy 0 ny jz 1 nz read 25 qz Jx JY 32 Jx 1 nx 3jy 1 ny 3jz 0 n2 eee ao mr E A ae 2 CN ts 2 Y q f 1 Steefel amp Yabusaki 46 OS3D GIMRT Manual or as read 23 qx jx jy jz jx 0 nx jy 1 ny jz 1 nz read 24 qay jx jy jz jx 1 nx jy 0 ny jz 1 nz for 2 2D problem i e where jz 1 The position in the driver routines os3d f and gimrt f to place an update of a transient velocity field is marked by c PLACE ROUTINE TO UPDATE TRANSIENT VELOCITY FIELD HERE 12 5 Adding a Temperature Field The appropriate place within the startup routine for GIMRT startup gi to place a steady state calculation of the temperature field is marked by the lines c PLACE ROUTINE TO INITIALIZE TEMPERATURE HERE BY READING IN C A TEMPERATURE FILE Cc Uncomment the call to the subroutine DENSITY c THIS WILL OVERWRITE ANYTHING SPECIFIED ABOVE Chkkkkkkkkkkkkkkkkkkkkkkkkkk kk c do jy 1 ny c do jx 1 nx c J jy 1 nx jx c call DENSITY j c end do c end do c kkkkkkkkkkkkkkkkkkkkkkkk k k
43. centrations but basis switching when carried out every Newton iteration does appear to reduce the number of times the Newton fails to converge particularly in the case of redox problems If full basis switching is used every Newton iteration it should not be necessary to solve for the logarithms of the concentrations and one can make use of standard relaxation techniques to reduce the size of the Newton step when the linear concentrations are employed The code carries out basis switching by first sorting the species concentrations in order of abundance and then systematically attempting to construct a linearly in dependent set from those concentrations see the subroutine switch The recommended procedure is to try to run the problem without basis switching since that will be the fastest method if the code con verges properly If numerous non convergences occur however we recommend trying basis switching The various input parameters required to turn the option on and off and to set the maximum allowable concentration correction limits are described in detail in the section on the input file threedin 5 7 Activity Coefficient Model At this time the user may choose to run the code either with unit activity coefficients or with an extended Debye Htickel formulation to calculate activity coefficients for the ionic species The Debye Hiickel formulation is given by AZ I log yi 01 77 1 Ba I The coefficient
44. ch Steefel and Lasaga 1994 Although the two codes do not yet have totally identical capabilities they can be run from the same input file allowing comparisons to be made between the two approaches in many cases The advantages and disadvantages of the two approaches are discussed more fully below but in general OS3D is designed for simulation of transient concentration fronts particularly under high Peclet number transport conditions because of its use of a total variation diminishing or TVD transport algorithm GIMRT is suited for simulating water rock al teration over long periods of time where the aqueous concentration field is at or close to a quasi stationary state and the numerical transport errors are less important Where water rock interaction occurs over ge ological periods of time GIMRT may be preferable to OS3D because of its ability to take larger time steps Both OS3D and GIMRT are based on a continuum representation of transport and reaction in porous media The codes solve the nonlinear set of partial differential equations describing coupled reaction and transport in a multicomponent system by using finite difference techniques The finite difference method although computationally more intensive than the simpler reaction path models allows for modeling of the entire range of both transient and steady state phenomena including advective dispersive and diffu sive transport or any combination of these A second important f
45. chemical condition but these will not remain fixed as transport and reactions act on the grid cells The geochem ical condition 1 is redundant since in the absence of information here the code will use condition 1 as the default but it is allowed nontheless This section is terminated when the code encounters a value for ndist 0 this must be there for a correct read JXXLO I grid point in X direction specifying lower left hand corner of a rectangular domain to be initialized with geochemical condition ndist JXXHI 1 grid point in X direction specifying upper right hand corner of a rectangular domain to be initialized with geochemical condition ndist JYYLO I grid point in Y direction specifying lower left hand corner of a rectangular domain to be initialized with geochemical condition ndist JYYHI grid point in Y direction specifying upper right hand corner of a rectangular domain to be initialized with geochemical condition ndist JZZLO I grid point in Z direction specifying lower left hand corner of a rectangular domain to be initialized with geochemical condition ndist JZZHI D grid point in Z direction specifying upper right hand corner of a rectangular domain to be initialized with geochemical condition ndist e Card 17 JBD I geochemical condition to be used at the boundaries of the 3D rectangular volume The code expects to read 6 numbers here one for each edge of the cube even for 2D 1D and 0D problems A neg
46. d the mineral volume fractions and mineral surface areas are updated If the reaction transport calculations are combined with a calculation of the flow and temperature field these would normally be updated as well at this time At the end of each time step an algorithm which com putes the second derivative of the solute concentrations with respect to time is used to determine whether the time step should be increased or decreased This algorithm provides a way of minimizing the time truncation error if that is an issue in a particular problem and also of controlling the size of the time step so as to maintain numerical stability The program structure for the time splitting operator splitting and sequential iteration method is similar differing primarily in the fact that the transport is carried out in advance of any calculations of the reactions eect orem nape CE oe tesen OS3D GIMRT Manual 19 Steefel amp Yabusaki Figure 3 The reaction calculations are then calculated individually at each point in space Since these reactions are not directly coupled via transport to the reactions in neighboring cells the calculations can be carried out independently The size of the system of equations to be solved at each grid point is then Ne by Ne and there is no need to assemble the global system of function residuals and the global Jacobian matrix over the entire spatial domain The sequential iteration method differs in
47. e Dirichlet boundary condition option is used usage described below The default is therefore no diffusive dispersive flux boundaries which is only overridden where the Dirichlet option is invoked 5 3 Choice of Primary and Secondary Species After setting up the problem physically one needs to decide on the chemistry to be included The method for specifying primary and secondary species is essentially that devised by Peter Lichtner for the code MPATH The codes OS3D and GIMRT require that one specify an initial choice of primary species which then determine the number of independent chemical components in the system This choice of pri mary species is not unique however as discussed in the section 2 2 By making use of the basis switching OS3D GIMRT Manual 23 Steefel amp Yabusaki option itis possible to allow the code to dynamically choose the set of primary species to be used in writ ing the reactions at each grid point see discussion below Unlike EQ3 EQ6 Wolery et al 1990 for instance the user must choose which secondary species are to be included in the simulation This can be an advantage for multicomponent reaction transport modeling since one would like the option of carrying out chemically simplified simulations However the user must take care that important species are not neglected Future versions of the code will have the capability of reading all of the potentially relevant species from the database The only
48. e added to the thermodynamic database e Card 9 SECONDARY SPECIES e NAMCX D name enclosed in single quotes of the species to be included yl e s marks end of secondary species read Steefel amp Yabusaki 32 OS3D GIMRT Manual 7 3 Minerals The reacting minerals to be included in the simulation are read in from card 10 with the user provided rate data An example of the format is given below 10 MINERAL K nreact k sati k sat2 k thresh k npre 11 k species dep mm 11 depend mm 11 k rate0 11 k Ea 11 k goethite 1 1 0 1 0 0 0 0 1 e 08 0 0 quartz 2 1 0 1 0 3 0 1 oh 0 3 6 31e 11 0 0 0 1 e 12 0 0 gibbsite 1 1 0 1 0 0 0 0 1 e 11 0 0 kaolinite 1 1 0 1 0 0 0 0 1 e 11 0 0 albite 3 1 0 1 0 0 0 1 ht 0 4 1 58e 10 0 0 1 oh 0 3 1 00e 10 0 0 0 1 58e 12 0 0 1 1 0 1 0 0 0 e Card 10 NAMK K NREACT K SAT1 K SAT2 K THRESH K NAMK K character string containing mineral name NREACT K number of parallel reactions involving the mineral The code will readnreact following lines SATI K order of reaction rate law n in Equation 32 SAT2 K order of reaction rate law M in Equation 32 THRESH K supersaturation log Q K necessary for precipitation e NPRE LL K DUMLABEL MM DEPEND MM LL K RATEO LL K EA LL K NPRE LL K number of species affecting the far from equilibrium dissolution rate The code will then read n
49. e appropriate thermodynamic data from for example the code SUPCRT The code also incorporates automatic basis switching rewriting the reactions in terms of the dominant species which can considerably improve the convergence behavior in some cases particularly redox or it can be run using a constant set of primary or basis species specified initially by the user The codes also allow for a range of mineral water reaction rate laws including parallel reactions e g pH dependent and pH independent rates and proton promoted and or ligand promoted dissolution At the present time only GIMRT can treat non isothermal conditions Neither code currently has a hydrodynamic flow calculation coupled directly to it At the present time therefore one must provide a velocity field calculated externally A fully coupled flow code with options Steefel amp Yabusaki 2 OS3D GIMRT Manual for porosity and permeability change as a function of time will be included in the next release 2 MATHEMATICAL FORMULATION 2 1 Conservation Equations The governing equation for the conservation of solute mass is given by Steefel 1992 Steefel and Lasaga 1994 0 z pf MmoC V e DV p Mp 0C upsMy 0C 1 1 2 Ntor where C is the molal concentration of a species in solution in units of moles per kg H20 My o is the mass fraction of H20 py is the fluid density u is the Darcian flux D is the combined dispersi
50. e codes OS3D and GIMRT the user has the option of either leav ing basis set always the same or of dynamically calculating a new basis set at each grid point in space so that the basis species chosen is the one with the largest concentration The basis switching makes for a more diagonally dominant Jacobian matrix and therefore may improve the convergence behavior of the algorithm 2 3 Reaction Rate Laws The codes OS3D and GIMRT use a kinetic rate law based on the assumption that attachment and detach ment of ions from mineral surfaces is the rate limiting step Steefel and Lasaga 1994 The principal computational advantage of this approach is that there are no additional concentrations of ions immedi ately adjacent to the mineral surfaces to be solved for as there are when diffusion to and from the surface is rate limiting Steefel 1992 Steefel and Lasaga 1994 Where surface attachment is the rate limiting step the concentrations next to the mineral s surface are just the concentrations of those ions in the bulk solution which we are already solving for If the mineral volume fractions are held constant for a single time step and updated only after convergence of the aqueous and surface species then we can solve for Miata raeme or Steefel amp Yabusaki 6 OS3D GIMRT Manual the minimum number of equations which is Ne x M where Ne is the number of independent chemical components and M is the number of g
51. eature of the codes is the use of a kinetic formulation for all mineral dissolution and precipitation reactions This allows one to avoid a priori as sumptions of local equilibrium although it is still possible to obtain local equilibrium behavior where the rates of reaction are much greater than the rates of mass transport In contrast aqueous and surface com plexation reactions including redox are assumed to be reversible at this time although it is planned that this requirement will be relaxed in future developments The package OS3D GIMRT includes a number of notable features OS3D can be run in OD D here refers to the spatial dimensions with OD being a reaction path or batch reaction calculation 1D 2D or 3D simply by changing the input file GIMRT possesses the same flexibility but is currently limited to two spatial dimensions At startup the codes allow one to specify a number of different zones within the domain of interest which have distinct initial solute concentrations mineral abundances and or grid spacing The code also includes a distribution of species calculation to specify the initial and boundary solute concentrations The code reads a thermodynamic data base based on the EQ3 EQ6 Wolery and others 1990 data base which covers temperatures between 0 and 300 C along the water saturation curve for a large number of mineral and aqueous species The user can simulate reactive transport at higher pressures by providing th
52. edox couple may not be important since the code will automatically write the reactions at any particular grid point in terms of the dominant redox couple i e that which is present in the largest concentration see section on Basis Switching below This may be particularly useful for redox problems which can show convergence problems this approach works best where basis switching is carried out every Newton iteration The user may choose not to use basis switching at all to use it at the beginning of every timestep or to use itevery Newton iteration The same proviso however applies about including ners p em m ou lt lt lt Steefel amp Yabusaki 24 OS3D GIMRT Manual the species which are used as basis species in the EQ3 EQ6 data base 5 5 Acid Base Reactions Acid base reactions require no special treatment other than in the initialization process where one nor mally fixes or calculates pH rather than an Ht concentration Once a pH has been either fixed or calculated in the initialization procedure however Ht can be treated like any other component Again no special considerations are needed in contrast to the statements by Yeh and Tripathi 1989 5 6 Basis Switching As discussed above both OS3D and GIMRT can rewrite the reactions in terms of the most concentrated set of linearly independent species This is not absolutely necessary if one is solving for the logarithms of the con
53. egrid a power law scheme proposed by Patankar 1980 is used This method slides between a fully centered form at Pegria lt 2 to an upstream weighted formulation at Pegria gt 10 where the upwind scheme is necessary for stability Using the power law scheme of Patankar 1980 the finite difference coefficients become 0 16 3 ag De 0 1 Es eet 0 52 e 0 15 aw Da o 1 LEI po nu 653 w taman A lt gt meee gs re gt ARA RITA OS3D GIMRT Manual 13 Steefel amp Yabusaki ap QGg ay 54 where the operator A B means that the greater of the values A or B 15 used equivalent to the FOR TRAN operator AMAXI A B and A refers to the absolute magnitude of the value A fuller descrip tion may be found in Patankar 1980 A fully implicit backward method is used to represent the time derivative Lapidus and Pinder 1982 The equations to be solved for each component at each grid point at the t At time level are then a UP t At UF t apUf t At agUF t At awUj t At RP At As 0 55 where At is the time step and t refers to the values at the present time level As discussed above the equations are written in terms of the concentrations of the primary species in the code GIMRT see Equa tion 32 and are expressed in terms of the total concentrations here in order to be concise 3 5 Transport Algorith
54. er of files which write out the values of the important quan tities e g pH reaction rate as a function of space at specified times during the course of the simulation Up to 19 files representing 19 different times may be created during any one simulation e g ph1 out ph19 out The files include e masterl out writes the concentration of the master variable specified in threedin If the species is H the code writes the solution pH Otherwise the logarithm of the concentration is written A OS3D GIMRT Manual 43 Steefel amp Yabusaki e concl out writes all of the linear concentrations e ratel out writes the reaction rates of the minerals in units of moles liter bulk volume sec7 The user may change the units to volume yr by uncommenting and commenting the appro priate lines in the graphics output files graphld f graph2d f and graph3d f e volumel out writes the volume of the minerals and the porosity e saturl out writes the saturation indices of the minerals log Q K The format of the output files may be easily modified by the user 11 PROBLEMS IN RUNNING THE CODE With a code as complicated as OS3D and GIMRT there are numerous things which could cause it to fail Potential users who want a seamless black box which solves their problems for them automatically are urged to consider a field other than reactive transport in porous media With some patience however it should be
55. erms of the user s choice of basis species We would also like to thank J E VanderKwaak P A Forsyth K T B MacQuarrie and E A Sudicky for making their WATSOLV sparse matrix solver available to us We are also pleased to acknowledge the support provided by the Geoscience Program of Basic Energy Sciences Dr William C Luth of the Department of Energy DE FG02 90 ER 14153 and to the Donors of the Petroleum Research Fund administered by the American Chemical Society PRF 22364 AC2 for the support of this research all to Antonio C Lasaga at Yale University In addition we would like to ac knowledge the support of the the Schweizerishes Nationalfonds in the form of a grant 21 30908 91 to Peter C Lichtner while at the Universit t Bern in Switzerland We are grateful for the financial support of Prof Dr Tjerk Peters at the Mineralogisch Petrographisches Institut at the Universit t Bern Finally we would like to thank Ashok Chilakapati of PNNL for his careful review of the manuscript OS3D GIMRT Manual 49 Steefel amp Yabusaki 14 REFERENCES Aagaard P and Helgeson H C 1982 Thermodynamic and kinetic constraints on reaction rates among minerals and aqueous solutions I Theoretical considerations Amer J Sci 282 237 285 Aris R 1965 Prolegomena to the rational analysis of systems of chemical reactions Arch Rational Mech and Anal 19 81 99 Bear J 1972 Dynamics of Fluids in Porous Media American Elsevier
56. es of the reactions rates are given by min dl ia SS vym mkm de uea 707 Ln a 65 where any dependence of the rate law on the activities of species far from equilibrium has been neglected and only first order reactions are considered M and n 1 For the case where the activities of various species e g H affect the far from equilibrium dissolution rate additional terms appear in the partial derivatives not given here The partial derivatives of the ionic strength J with respect to the primary species can be written as Lichtner 1992 91 za 2 D227 Cr 2Cy 907 11 235 22 12 vadla7 d diny dl and the derivatives of the activity coefficients with respect to the 10 16 strength are dlog yk _ Az 3 lt lt b 67 211 2 1 B ak Pe If we define f as making up the elements of the vector b when taken over all of the grid points in the system and the 9 f OC as forming the elements of the global Jacobian matrix A then Equation 59 can be written in the familiar form A 0C b 68 Once the dC are computed they are used to update the concentrations of the primary species 0 C3 400 4 1 Ne x M 69 This completes a single Newton iteration The procedure is repeated until convergence of the function residuals the f s are reduced to a tolerance set in the code i e fixo 70 Without basis switching it may be necessary to use the logarithms of the concentrations since
57. espond by suggesting that the user fix an interior node using Card 16 e Card 16 NDIST JIXXLO D JXXHI D JY YLO D JY YHI D JZZLO D JZZHI D NDIST I geochemical condition number according to the order given in Card 15 to be distributed in space If the geochemical number is negative see second example above the code will treat this domain as fixed throughout the simulation If the geochemical condition is positive then the code will initialize the domain to the appropriate geochemical condition but these will not remain fixed as transport and reactions act on the grid cells The geochemical condition 1 is redundant since in the absence of information here the code will use condition 1 as the default but it is allowed nontheless This section is terminated when the code encounters a value for ndist 0 this must be there for a correct read JXXLO D grid point in X direction specifying lower left hand corner of a rectangular do main to be initialized with geochemical condition ndist JXXHI grid point in X direction specifying upper right hand corner of a rectangular do main to be initialized with geochemical condition ndist JYYLO D grid point in Y direction specifying lower left hand corner of a rectangular do main to be initialized with geochemical condition ndist JYYHI D grid point in Y direction specifying upper right hand corner of a rectangular do main to be initialized with geochemical condition ndist JZZLO g
58. fficients 26 INSTALLATION st pra Ed aya Rie E EA ARA RRA A a 26 RUNNING GIMRT and OS3D ia o A ws A a a 28 71 Control P rameterS oros a ae 28 7 2 Specifying Primary and Secondary Species e 30 13 Minerals Sorama e eusen rs Deed Gi Sed ee eet ee Ses E S 32 T4 A te Boas ee as Bae E OE 33 7 5 Discretization 4c ori Ps E De eta ee Se a SS ee SESS 33 7 6 Geochemical Conditions pa bao ee e Sie caw wig Bi ge 34 7 7 Spatial Distribution of Geochemical Conditions 36 7 8 Reaction Path Switch and Plotting Information 39 7 9 Transport Patameters cit a a ede A SAA en ee 39 PARAMETERS IN FORTRAN FILES 00 4 OUTPUT IN THE FILE THREEDOUT 5 AAA 42 GRAPHIESQUIPUT aeie Ta E IIA a ee A ORE E a 42 PROBLEMS IN RUNNING THE CODE 2 0 0 43 CUSTOMIZING THE CODE 300540 A A Sa 44 12 1 Array Dimensions zarete e E a a eg oe 44 12 2 Modifying the Data Base a AAA A a te ASA 44 12 3 Adding or Modifying Kinetic Rate Laws oscs ms ss 45 12 4 Providing and or Updating a Velocity Field 45 12 5 Adding a Temperature Field os 2 A Bee Bea ee 46 OS3D GIMRT Manual ili Steefel amp Yabusaki 13 ACKNOWLEDGMENTS
59. for the treatment of its activity coefficient OS3D and GIMRT assume that the activity of H20 is given by its mole fraction in the solution 5 4 Redox Reactions There are many possibilities for choices of primary and secondary species to represent redox reactions Both OS3D and GIMRT write all redox reactions as whole cell reactions i e electrons are not used as either primary or secondary species It is easy to show that the electrons are completely unnecessary if the solutions are electrically neutral and the inclusion of electrons requires that one include a species the mass of which should everywhere be zero With the whole cell reactions no special provision for redox species and reactions is necessary in contrast to statements by Yeh and Tripathi 1989 although the inclusion of redox reactions implies that we have included an additional mass balance equation representing the conservation of exchangeable charge Since Og is used as the basis species in redox reactions in the EQ3 EQ6 data base it must be included in the list of gases in the input file if redox reactions are to be considered A gas however cannot be used as a primary species in OS3D and GIMRT since gases are not included in the mass balances One could use O2z aq however or any number of other species which form a redox couple for example Fet and Fe 3 or HS and SO 2 When using the basis switching option the initial choice of primary species to represent a r
60. he destruction of one mole of the secondary species Equation 8 implies that the rate of production of a primary component y due to homogeneous reactions can be written in terms of the sum of the total rates of production of the secondary species Kirkner and Reeves 1988 Nz i DIAN 10 1 1 where r are the reaction rates of the secondary species Equation 10 suggests that one can think of a mineral dissolving for example as producing only primary species which then equilibrate instantly with the secondary species in the system Using Equation 10 the rates of the equilibrium reactions can be eliminated Lichtner 1985 Kirkner and Reeves 1988 2 05 Mmo c EA 14 X Ve ups Mazo c Da vX DV pr Mmo c SRE Ue ls Ne 11 S 7 ee OS Steefel amp Yabusaki 4 OS3D GIMRT Manual Note that only the term Roa remains on the right hand side of Equation 11 because we have assumed that they are the only kinetic reactions 2 2 Example of a Total Concentration If a total concentration U is defined Reed 1982 Lichtner 1985 Kirkner and Reeves 1988 Nz U C Y 17 X 12 1 1 then the governing differential equations can be written in terms of the total concentrations in the case where only aqueous and therefore mobile species are involved Kirkner and Reeves 1988 min J p MH003 V e up Mp 0U0 DV
61. hemical reaction Model development and application to water rock interaction Ph D thesis New Haven Connecticut Yale University 234 p Steefel C I and Lasaga A C 1990 Evolution of dissolution patterns Permeability change due to cou pled flow and reaction In Chemical Modeling in Aqueous Systems II ed D C Melchior and R L Bassett ACS Symp Ser No 416 212 225 Steefel C I and Lasaga A C 1992 Putting transport into water rock interaction models Geology 20 680 684 Steefel C I and Lasaga A C 1994 A coupled model for transport of multiple chemical species and ki netic precipitation dissolution reactions with application to reactive flow in single phase hydrother mal systems American Journal of Science 294 529 592 Steefel C I and Lichtner P C 1994 Diffusion and reaction in rock matrix bordering a hyperalkaline fluid filled fracture Geochim Cosmochim Acta 58 3595 3612 Steefel C I and Van Cappellen P 1990 A new kinetic approach to modeling water rock interaction The role of nucleation precursors and Ostwald ripening Geochim Cosmochim Acta 54 2657 2677 Strang G 1968 On the construction and comparison of difference schemes SIAM J Numer Anal 5 506 517 2 om A A w CA SO Meee OS3D GIMRT Manual 51 Steefel amp Yabusaki Valocchi Albert J and Malmstead Michael 1992 Accuracy of operator splitting for advection dispersion reaction problems Water Resources Res bf 28
62. hen read npre character strings dumlabe1 and exponential dependences depend DUMLABEL MM character string containing species name affecting far from equilibrium disso lution rate DEPEND MM LL K exponential dependence of the far from equilibrium rate on species given in preceding dumlabel RATEO LL K far from equilibrium mineral dissolution rate constant at 25 C in units of moles m s 1 Can be used as the rate constant at temperature if activation energy is set 0 EA LL K activation energy in units of kcal mole end of mineral read woe Steefel amp Yabusaki 54 OS3D GIMRT Manual e Card 11 NAMG D name enclosed in single quotes of the gas to be included 1 o s marks end of read e Card 12 NZONEX NVX D DXX NZONEX numter of zones in X direction with different grid spacing NVX I number of grid cells using the grid spacing DXX D grid spacing e Card 13 NZONEY NVY I DYY D NZONEY number of zones in Y direction with different grid spacing NVY D number of grid cells using the grid spacing DYY D grid spacing in Y direction e Card 14 NZONEZ NVZ D DZZ NZONEZ number of zones in z direction with different grid spacing NVZ D number of grid cells using the grid spacing DZZ I grid spacing e Card 15 NCHEM number of geochemically distinct initial and or boundary conditions to be included In what follows the subscript K refers to the geochemical conditi
63. in the system specified in threedin gt e mrct 1 plus the number of reacting minerals in the system specified in threedin e mspec 1 plus the number of secondary species in the system specified in threedin e mgas plus the number of gases in the system specified in threedin and e mx my mz the number of grid points in the X Y and Z directions respectively used to solve the for the solute concentrations specified in threedin If the file params inc or params os inc are modified then all of the subroutines which include this file are recompiled when Makefile is run again When the appropriate changes have been made the source code can be compiled and an executable image can be created This is done by giving the command make which runs the default file Makefile So as to avoid the necessity of creating multiple copies of the executable image one normally creates a symbolic link to the user s bin directory This allows the code to be run from any directory chosen by the user assuming the input file threedin is located there This can be done by moving to the user s bin directory with the command cd bin amm a a rn E O 2 Vea a ae ae E Steefel amp Yabusaki 28 OS3D GIMRT Manual where indicates the user s home directory and then typing In s path gimrt gimrt followed by rehash to update the hash table so that the symbolic link is recognized The designation path refers to the lo
64. inematic dispersion coefficient is written here as uz D arlul ay or ea 6 where az and ar are the longitudinal and transverse dispersivities respectively Bear 1979 Marsily 1986 and u is the magnitude of the velocity In both OS3D and GIMRT off diagonal dispersion co efficients i e Diz where 1 and 7 are coordinate directions are not included at the present time eee eet lt LE tee lt lt lt n ee 7 gt OS3D GIMRT Manual 3 Steefel amp Yabusaki The heterogeneous reaction term R7 can be written as the sum of all the individual mineral water re actions which affect the concentration of the zth species Nm REA a e UmTm 7 m 1 where rm is the rate of precipitation or dissolution of mineral m per unit volume rock V m is the number of moles of mineral m or stoichiometric coefficient if the reaction is written in terms of the dissolution of one mole of the mineral and Nm is the number of minerals present in the rock Following convention we take rm as positive for precipitation and negative for dissolution Equation 1 provides a general formulation for the conservation of solute mass which makes no assump tions of chemical equilibrium If we assume that the various aqueous species are in chemical equilib rium however it is possible to reduce the number of independent concentrations that is the number that actually need to be solved for Mathematically th
65. is 25 C or tomastertemp data if another temperature is given e Card 3 TSTEP DELT TTOL JACOB BPEG TSTEP maximum time step in years allowed DELT initial time step in years TTOL parameter controlling the tolerance on increasing the size of the time step Making ttol smaller decreases the rate at which the time step is increased for the purpose of controlling either the time truncation error or for numerical stability IJACOB flag controlling update of Jacobian 0 update Jacobian every Newton iteration 1 whendelt 1 2 tstep update Jacobian every 3 Newton iterations and when delt tstep update Jacobian every 6 iterations BPEG maximum change in molality in any one Newton iteration Reducing bpeg increases the under relaxation of the Newton step and may sometimes prevent nonconvergence although the total number of iterations required may actually increase e Card 4 JPOR CONSTANTPOR VDISSMAX VPPTMAX JPOR flag controlling update of porosity and flow O no update of porosity 1 porosity is updated CONSTANTPOR porosity of the medium for the case where the porosity is not updated jpor 0 If jpor 0 the porosity is taken from constantpor rather than from the sum of the volume fractions of the minerals in the system If jpor 1 porosity updated constantpor is ignored VDISSMAX maximum rate of decrease of a mineral s volume fraction units of yr Can be used to limit the amount of mineral di
66. is means that in a system containing N o aqueous species the number of independent chemical components in the system N is reduced from the total num ber of species by the Ny linearly independent chemical reactions between them for further discussion see Hooyman 1961 Aris 1965 Bowen 1968 Van Zeggeren and Storey 1970 Reed 1982 Lichtner 1985 Kirkner and Reeves 1988 This leads to a natural partitioning of the system into Ne primary or ba sis species designated here as C and the N secondary species referred to as X Reed 1982 Lichtner 1985 Kirkner and Reeves 1988 The equilibrium chemical reactions between the primary and secondary species take the form Ne 4 So Ay 1 Nz 8 j 1 where the A and the A are the chemical formulas of the primary and secondary species respectively and is the number of moles of primary species 7 in one mole of secondary species 2 It should be noted here that the partitioning between the primary and secondary species is not unique that is we can write the chemical reactions in more than one way The reversible reactions provide an algebraic link between the primary and secondary species via the law of mass action for each reaction Ne X Key mo 2 1 snag Nz 9 j 1 where y and y are the activity coefficients for the primary and secondary species respectively and the K are the equilibrium constants of the reaction given in Equation 8 written here as t
67. ity in a particular coordinate di rection or to give 3 separate filenames where the code can read the X Y and Z components of the velocity A AN ms P OA lt x ES tt A Steefel amp Yabusaki 40 OS3D GIMRT Manual field Note that the code will attempt to read 3 character strings for the velocity files even if the parameter ivelocity is set to 0 In the same way the code will read the components of the constant velocity field even if this information is not going to be used After this information is read the code reads the value of the molecular diffusion coefficient to be used A sample input follows 20 VELOCITY FIELD OPTIONS 0 CONSTANT VELOCITIES 1 FILE READ in the order QX OY 02 0 0 5 0 0 0 0 custom32x16 vx custom32x16 vy custom32x16 vz 21 HYDRODYNAMIC DISPERSION 0 000000001 0 000000001 22 MOLECULAR DIFFUSION idiffus dcoeff porpow 1 1 d 09 1 0 e Card 20 IVELOCITY flag indicating whether to use a constant velocity field or to use the information found in the files specified below 0 assume constant velocities given in gxinit qyinit and qzinit below 1 use velocity components found in the files listed below e QXINIT QYINIT QZINIT X Y and Z components of the constant velocity field Code will attempt to read these numbers even if ivelocity 1 but the information will not be used e VXFILE VYFILE VZFILE character strings containing file names with
68. kkx k If one wishes to include a transient calculation of the temperature field the position in gimrt f inside the time loop is indicated by the lines c PLACE ROUTINE TO UPDATE TRANSIENT TEMPERATURE FIELD HERE 6 Uncomment the call to the subroutine DENSITY cChkkkkkkkkkkkk kkkk kkxkkkk kkkkkk c do jy 1 ny do jx 1 nx c J jy 1 nx jx c call DENSITY j c end do c end do c kkkkkkkkkkkkkkkk kkxkkkkkk k kkkkx k The temperature field is block centered unlike the face centered temperatures so the temperature will have dimensions of 1 nx 1 ny 1 nz gt gt gt ae OS3D GIMRT Manual 47 Steefel amp Yabusaki Table 1 Sample data base read by GIMRT and OS3D for a single temperature of 25 C temp 1 25 h20 3 0 0 0 18 01528 al 3 9 0 3 0 26 98154 9 0 1 0 1 00794 hco3 4 0 1 0 61 01714 sio2 aq 3 0 0 0 60 08430 null 0 0 0 0 al oh 2 3 2 0 h 1 0 al 3 2 0 h2o 10 0991 4 00 1 00 60 99622 al oh 3 aq 3 3 0 h 1 0 al 3 3 0 h20 16 1577 3 00 0 00 78 00356 al oh 4 3 4 0 h 1 0 al 3 4 0 h20 22 1477 4 00 1 00 95 01090 aloh 2 3 1 0 h 1 0 al 3 1 0 h20 5 0114 4 50 2 00 43 98888 602 3 1 0 h20 1 0 ht 1 0 hco3 6 3447 3 00 0 00 44 00980 co3 2 2 1 0 h 1 0 heo3 10 3288 5 00 2 00 60 00920 null
69. lly changing surface area of the mineral e NAMK D VOLIN K AREAIN K NAMK D character string containing mineral name VOLIN K initial volume fraction of mineral I in the Kth zone AREAIN K initial reactive surface area of mineral I in the Kth zone Units of M in ral 3 m rock 7 7 Spatial Distribution of Geochemical Conditions Information on the spatial distribution of geochemical conditions is read in Cards 16 and 17 The code assumes that geochemical condition 1 is the default and will initialize the entire interior domain with these concentrations and constraints if no other internal zones are specified in Card 16 The internal domains are assumed to be orthogonal blocks which may have dimensions of zero width In the example which follows one internal zone is specified which begins at JX 1 and ends at JX 3 inclusive This example could apply either to a 1D problem where we automatically assume NY 1 or it could be a 2D problem where we are initializing a strip 3 grid cells along JY 1 The code will read up to 100 internal zones i e the number of heterogeneities and will terminate only when it finds a geochemical condition equal to 0 the geochemical condition 0 in this respect is similar to used above to terminate character string reads This section can also be used to specify the kind of boundary condition and or initial condition desired in addition to the geochemical condition to be used In
70. lt species name gt lt molar volume gt lt of component species in reaction gt lt stoichiometric coeff gt lt component name gt lt logK s N values gt lt molecular weight gt lt comment gt The sign of the log K must conform to the reaction written with the aqueous complex mineral or gas on the left hand side that is as the destruction of one mole of the aqueous complex mineral or gas A sample data base is shown in Table 12 5 The data base delivered with the code is a transformed version of the EQ3 6 data base created by Peter Lichtner at Southwest Research Institute This data base contains over 1500 aqueous species and minerals Log K values are stored at temperatures of 0 25 60 100 150 200 250 and 300 C and interpolated according to the expression oe Oe ln ha ost 79 logK T2 T where T denotes the absolute temperature in degrees Kelvin Pressure lies along the steam saturation curve Other data bases may be generated using the computer code SUPCRT Johnson et al 1992 for any desired temperature range and pressure 12 3 Adding or Modifying Kinetic Rate Laws The kinetic rate laws presently implemented in GIMRT and OS3D all include some kind of saturation dependence of the rate This ensures that the rates go to zero when equilibrium is reached Currently the codes assume the rate law prevails during precipitation although one can specify obtain some differ ences in the dissolution and
71. m in Operator Splitting Approach 3 5 1 Advection When solutions require accurate resolution of concentration fronts in advection dominant problems high resolution shock capturing techniques are often employed These solution techniques are typically based on explicit numerical schemes that are incompatible with a globally implicit solution of coupled trans port and reactions e g GIMRT Decoupling i e time splitting the transport algorithm from the geo chemical reaction algorithm relaxes the requirement for an implicit transport method allowing the use of high resolution shock capturing techniques based on explicit numerical schemes The principal trade off for the additional accuracy is a strict limit on the ratio of time step to grid spacing a limit that is in herent in explicit transport methods Typically this limit becomes most restrictive under steady state or near steady state conditions Note that steady state does not mean static processes In physical systems steady state occurs when dynamic processes are balanced or fully compensating In OS3D an explicit third order total variation diminishing TVD numerical scheme Gupta et al 1991 is used to model non dispersive advective transport The technique is based on calculating an estimated concentration at the face of a grid cell that is dependent on the flow direction In the case of flow in the positive x direction i e gt 0 the estimated concentration would be
72. meters which may be used are described in detail in the WATSOLV manual but they are reviewed briefly here The data statement at the top of the files gimrt f and newton f take the form watsolv input parameters data level 1 idetail 0 rtwotol 1 e 25 rmaxtol 1 e 25 amp smaxtol 1 e 12 nitmax 100 solve_flag false north 2 The parameters have the following meaning level factorization level number of grid connections idetail solver performance detail level O no information 1 summary information 2 full information rtwotol residual Euclidean two norm convergence error criteria rmaxtol residual infinity maximum norm convergence error criteria smaxtol solution update infinity maximum norm scale criteria nitmax maximum number of iterations to perform solve_flag logical flag to indicate acceleration method true GMRES acceleration false CGSTAB acceleration north number of basis vectors to retain in GMRES Other input parameters for GIMRT are also listed in data statements e startup parameters data nend 5000000 newton 50 corrmax 5 0 atol 1 e 09 rtol 1 e 07 These parameters are defined as nend total number of allowable timesteps newton maximum number of Newton iterations corrmax maximum correction in log space for Newton update atol absolute tolerance for convergence on residuals mol m ttol relative tolerance for convergence on residuals mol m For OS3D there are some additio
73. mount of physical memory gt 16 megabytes 12 2 Modifying the Data Base The data base read by GIMRT and OS3D is divided into four blocks The first block is a list of basis species not included on the left hand side of a chemical reaction with the 00 radius parameter used in the Debye Hiickel equation the species valence and the molecular weight of the species The second block contains chemical reactions and thermodynamic log K data for aqueous species the third block mineral data and the fourth block data for gases The various blocks are separated by lines beginning with the string null The first line of the data base contains the temperature for each log K entry with the total number of values First line temp NT 15 Ty Format for basis species lt component name gt lt ion size gt lt charge gt lt molecular weight gt Format for aqueous complexes lt species name gt lt of component species in reaction gt lt stoichiometric coeff gt lt component name gt lt logK s N values gt lt ion size gt lt charge gt lt molecular weight gt lt comment gt Format for minerals lt species name gt lt molar volume gt lt of component species in reaction gt lt stoichiometric coeff gt lt component name gt lt logK s N values gt lt molecular weight gt lt comment gt OS3D GIMRT Manual 45 Steefel amp Yabusaki Format for gases
74. n The next part lists the primary and secondary species the minerals and the gases chosen for the run The code then writes out the stoichiometric coefficient ma trix and the log Ks for the secondary species minerals and gases The user would do well to examine this matrix since it represents in tableau form the reactions written in terms of the user provided pri mary species The code has also transformed the log Ks in the EQ3 database into the appropriate log K for the reaction as written given the user provided primary species Following this is information on the mineral rate laws basically summarizing what the user provided in Card 10 of threedin Next is the information on the discretization and the geochemical conditions provided at this point just regurgitation of the user provided input Below this section is information on the surface complexation including the calculated total concentrations of surface hydroxyls on the various minerals in the system Then the code proceeds to the actual initialization routine and provides full speciation given the constraints provided by the user This section includes information on the ionic strength the pH the charge balance and the concentration of the individual primary and secondary species after the convergence of the initialization routine At the end is a summary of the transport parameters which will be used in the run 10 GRAPHICS OUTPUT The subroutine graphics f creates a numb
75. n reduce the value of t tol in the input file threedin This parameter controls the rate at which the time step is increased based on the time truncation error In addition one can specify a smaller value for the maximum timestep tstep If the code blows up during the first few time steps then try reducing the size of the initial time step Once the system relaxes the time step will be gradually increased Non convergence during the time stepping rou tine is generally more likely as the reaction rates are increased toward the local equilibrium limit Large reaction rate constants may require the use of smaller time steps If necessary one may also try a dif ferent value for the parameter bpeg if basis switching is used Bpeg may be changed in the input file threedin 12 CUSTOMIZING THE CODE There are a number of changes which can be made either to the code or to the data bases to customize GIMRT and OS3D for a particular problem 12 1 Array Dimensions As mentioned above it may be necessary or advisable to change the dimensions used in the code which are specified in the file params inc and paramsos inc If this file is changed GIMRT and OS3D must be recompiled with the make command The code checks to make sure that the array dimensions are large enough but it will not of course scale the dimensions back for smaller problems The dimension of the arrays may be an issue if GIMRT in particular is being run on a machine without a fair a
76. nal parameters set in the driver routine os3d the most important of which are logical flags choosing whether the code should be run with sequential iteration or straight operator splitting i e sequential coupling These include cee gape pee song qa a he ga a Cee ENS MESS Steefel amp Yabusaki 42 OS3D GIMRT Manual data nend 5000000 picard false picardtol 1 e 10 omeg 0 50 where nend total number of allowable timesteps picard logical flag indicating whether to use sequential iteration or coupling false use classic operator splitting or sequential coupling true use sequential iteration picardtol tolerance mol kg on convergence of sequential iteration omeg damping of sequential iteration to improve convergence We recommend the use of omeg equal to 0 5 or the sequential iteration method does not converge partic ularly well In the file newton 1 the user will find parameters similar to those listed under gimrt above controlling the number of allowable Newton iterations and the convergence tolerances 9 OUTPUT IN THE FILE THREEDOUT Upon completion of the initialization procedure for either OS3D or GIMRT information about the ini tialization the stoichiometric reaction matrices i e how the reactions are actually written the log Ks and a host of other parameters is written out to the file threedout The first part of the file summarizes the various control parameters used in the ru
77. nd within any one reaction mechanism multiple dependences on species far from equilibrium All reactions however include a saturation state dependence which guarantees that the rate will go to zero when the mineral reaches equilibrium with the solution 7 4 Gases The format for inputting gases which are not used in the mass balances is the same as the format for secondary species and are read from Card 11 e Card 11 NAMG D name enclosed in single quotes of the gas to be included ely marks end of read 7 5 Discretization The next section of the input file threedin Cards 12 through 14 includes user provided information about the discretization to be used in the problem All 3 dimensions must be specified even for 2D 1D or OD problems since the code will expect to read three lines here An example of a 1D problem using 32 grid cells in the X direction with a grid spacing of 1 0 meter follows Note that since the number of zones in the Y and Z directions are set to O the information about the number of grid cells in the Y and Z directions is ignored This problem would be converted to a 2D problem simply by changing nzoney to 1 thus yielding a 32 by 16 rectangular regime with grid cells measuring 1 meter by 1 meter If the code is to be run in OD i e a reaction path calculation then nzonex nzoney and nzonez are all set to 0 This gives nx ny and nz all equal to 1 12 nzonex nvx 1 dxx i i 1 nzonex 1 32 1 0
78. non As such the formulation neglects many fea tures that have been observed in real systems Principally the user should be aware that unlike molecular diffusion dispersion is length and time scale dependent because of the inhomogeneity of the physical media As plumes start to spread they are exposed to larger and larger scales of spatially variable ma terial properties Consequently it is impossible to define dispersion parameters that are effective for all time Moreover the magnitude of the dispersion parameter is dependent upon the grid resolution larger for large grid cells and smaller for small grid cells 3 6 Solving the Equations with Newton s Method As an example consider a problem involving one dimensional multicomponent reactive transport Equa tion 55 for any component 7 at the grid point P is a nonlinear function of the concentrations of the pri mary species at the the P grid cell in the case of OS3D and of the P E and W nodal points in the case of GIMRT An iterative scheme therefore is required to solve the set of nonlinear equations in either case Both codes solve the nonlinear equations with Newton s method which makes use of a Taylor se ries expansion to linearize the equations Using Newton s method the linearized set of equations for the 1D global implicit case become Ne af 3 Ne 0 P E Ne are 350 J fP 2o 0210 gt cE 2 oro HE 59 y 1 1 where fj is the discretized diffe
79. ns these represent the conditions which are geochemically distinct These include information on the constraints to be used in initializing the concentrations and initial mineral volume per centages and surface areas The list includes information on all of the primary species and minerals spec ified above in the input file threedin the code checks to make sure the lists match the list in Card 8 In this case the code will not look for an asterisk to terminate the read since it is expecting ncomp lines with information on the primary species and nrct lines on the minerals This part of the input file takes the form 15 nchem number of distinct geochemical conditions i 1 ncomp k 1 nchem ulab i itype i k guess i k ctot i k ncon i k OS3D GIMRT Manual 60 2 tracer h fe 3 al 3 sio2 aq nat 1 hco3 gt feoh goethite quartz gibbsite kaolinite albite 60 2 tracer h fet 3 al 3 102 80 nat 1 hco3 gt feoh goethite quartz gibbsite kaolinite albite 001 25 00 00 65 CTO COO URN WW WW 001 25 00 00 65 OOOO OM RN PP A 235 e 10 0e 10 0e 9 OE 15 0e 20 0 4 5be 1 5e 1 0e 5 e 6 100 0 10 0 10 0 10 0 10 0 e 05 001 0e 5
80. ogeneous reactions and co refers to the transported and partly reacted intermediate concentrations This is fol lowed by _ Um Note that the time step for reaction step is only 1 2 of the full timestep used in the transport step This results in a reaction term which is centered in time Walter et al 1994 The constant reaction term used in Equation 38 is calculated at the end of each reaction step from 2 opt 2 Re At 2 40 Since the term Ry used in Equation 38 is only approximate iteration is required between Equation 38 and Equation 39 The sequential iteration method often does not change the result obtained using classic time splitting assuming the Courant number is about 0 2 or less The user of the code may verify this for themselves with their own application by running OS3D in both modes procedure described below 3 3 Integrated Finite Difference Formulation The set of partial differential equations represented by Equation 32 is discretized using the integrated finite difference method Patankar 1980 Marsily 1986 As an example of how the discretization pro ceeds consider a set of partial differential equations for the total soluble concentration where the only important transport is via diffusion in one dimension and for simplicity we restrict ourselves to steady state conditions and to a system with constant density The differential equation is then OU D 2 RF 0 Oz s O
81. on diffusion tensor R in units of moles per unit volume rock per unit time is the total reaction rate of species 2 in solution and Nzoz is the total number of aqueous species Note that the porosity is imbedded in both the flux terms and the reaction terms as described below The reaction term R is divided into dissolution precipitation heterogeneous reactions R aqueous homogeneous reactions R and adsorption reactions R 4 such that R RP RY REY 2 The reaction rate is written here in terms of a unit volume of rock but could be written equivalently per unit volume of fluid if it were then multiplied by the porosity The Darcian flux is related to the true velocity of the fluid v by u v 3 D is a second order tensor which is defined as the sum of the mechanical or kinematic dispersion coeffi cient D and the molecular diffusion coefficient Do in water divided by the formation resistivity factor F x Do D D 4 F is defined as the ratio of the resistivity of the saturated porous medium over the resistivity of the pore solution alone Marsily 1986 The resistivity factor F may be defined in numerous ways see for ex ample Dullien 1979 but here a definition based on Archie s Law is used which gives the formation factor as F 5 where m is the cementation exponent Dullien 1979 reports values of m which range between about 1 3 and 2 5 The k
82. on number while I refers to the primary species e NAMC ITYPE 1 K GUESS I K CTOT 1 K NCON I K NAMC D name of primary species which must match that specified in Card 8 ITYPE K constraint to be used in initialization 1 mass balance using total concentration ctot 2 charge balance using the named species 3 mineral specified by ncon equilibrium constraint 4 pressure of gas bars fixed by ctot j i 5 species 15 a surface complex developed on the mineral specified by ncon In this special case guess and ctot taken on different meanings Guess in the case of a surface complex gives the site density on the mineral in units of moles surface hydroxyls per m mineral Ctot gives the intrinsic surface area of the mineral to be used in computing an initial total number of surface hydroxyl sites per mass solution units of m g7 7 activity of species is fixed by guess j i 8 concentration of species is fixed by guess j i GUESS I K initial guess for species concentration In the case of options 7 or 8 guess actually fixes either the concentration or the activity of the species In the case of option 5 pri mary species is a surface complex guess is the moles surface hydroxyls per m min eral CTOT 1 K total concentration of species if itype 1 If itype 4 gas constraint ctot is used to specify the partial pressure of the gas e g Oo while guess is used as a guess of the concent
83. or The global implicit approach of course eliminates this error but at the cost of adding in transport errors associated with the spatial dis cretization An alternative is to iterate sequentially between the transport and reaction in order to obtain a solution free from operator splitting error The sequential iteration method offers a way of doing this while still avoiding the construction and solution of a large global matrix as in the global implicit ap proach Our experience however is that the sequential iteration scheme offers little improvement over the Strang operating splitting procedure and that it involves considerably more CPU time At the present Steefel amp Yabusaki 10 OS3D GIMRT Manual time therefore we recommend using the Strang operator splitting scheme rather than sequential itera tion We have incorporated the sequential iteration approach nontheless and describe it below The sequential iteration approach is clearly explained by Walter et al 1994 In this case one solves the physical transport equations with a constant reaction term based on species concentrations from the previous iteration Each transport plus constant reaction term calculation is then followed by a reaction step using the transported and partly reacted concentrations The algorithm then becomes CC a LCi Re k 1 No No 38 where Ry is the total reaction rate for the individual species i e both heterogeneous and hom
84. ortions of this document may be illegible in electronic image products Images are produced from the best available original document Steefel amp Yabusaki 2 OS3D GIMRT Manual Copyright 1995 1996 by Battelle Memorial Institute All Rights Reserved OS3D GIMRT IS PROVIDED AS IS AND WITHOUT ANY WARRANTY EXPRESS OR IMPLIED THE USER ASSUMES ALL RISKS OF USING OS3D GIMRT THERE IS NO CLAIM OF THE MER CHANTABILITY OR FITNESS FOR A PARTICULAR PURPOSE YOU MAY MODIFY THE SOURCE CODE FOR YOUR OWN USE BUT YOU MAY NOT DIS TRIBUTE EITHER THE ORIGINAL OR THE MODIFIED CODE TO USERS AT ANY SITES OTHER THAN YOUR OWN OS3D GIMRT Manual i Steefel amp Yabusaki Contents LIST OE CONTENTS usas tea e a Sd Ee Bee a i E INTRODUCTION ans a id sd id GeO cogs Be 2 MATHEMATICAL FORMULATION eke ee ee e 2 2 1 Conservation Equations 2 2 2 Example of a Total Concentration il a i oe BSS es 4 2 3 Reaction Rate Laws isse weg hog ee dk Re a 5 2 4 Example involving Albite Dissolution lt 7 3 NUMERICAL FORMULATION 4 0 se ental www rd a 7 3 1 Global Implicit Approach 8 3 2 Operator Splitting Sequential Iteration Approach 9 3 3 Integrated Finite Difference Formulation no e ee SA Oe A 10 3 4 Transport Algorithm in Global Implicit App
85. p Mx 0U 9 1 Ne 13 As pointed out by Reed 1982 and Lichtner 1985 the total concentrations can usually be interpreted in a straightforward fashion as the total elemental concentrations e g total aluminum in solution but in the case of H and redox species the total concentration has no simple physical meaning and the to tal concentrations may take on negative values Where the system includes sorbed species i e surface complexes which are immobile then it is a necessary to define both a mobile total concentration Nz Z gt ARRE 14 wl and an immobile total concentration made up of the surface complexes Nz immobile _ grimmobile immobile U Ci So vX x 15 1 1 The governing differential equation then becomes 2 pp My o Ujrobile Ges 16 Ve ups My 0 Ure lt DV p MyoU3rvte 4 Run Melo No The concept of a total concentration will be clarified with a few examples Consider a system in which the aqueous species consist of H20 Ht OH AIF and Al OH Assuming that these species are all in equilibrium then the number of unknown concentrations i e Ne is reduced from the total number of species concentrations by the number of independent equilibrium reactions Nz between them Alt 4H O 4 0 4H 17 Ht OH H0 18 Choosing H20 Ht and Al as the primary species and writing the chemical reactions in matrix form as in Equation 8 tha
86. path 0 for batch 1 19 PLOTTING INFO AT SPECIFIC POINTS IN SPACE for reaction path calculations and breakthrough curves FIRST READ THE SPECIES NAMES TO WRITTEN TO OUTPUT IF NPLOT gt 0 THEN READ INTERVAL JXPLOT JYPLOT JZPLOT h heo3 kt 2 1 1 1 e Card 18 IPATH parameter only used when nxyz I i e no spatial discretization and controlling whether the mineral volume fractions are updated When ipath 1 and nxyz 1 the code runs in reaction path mode in order to model the movement of a fluid packet through the rock leaving reaction products behind When ipath 0 the code runs in true batch mode and the reaction products remain in the closed system e Card 19 PLOTDUMMY I character string for species to be plotted at a specific point in space i marks end of read e INTERVAL JXPLOT JYPLOT JZPLOT interval at which time dependent data the species listed above in plotdummy will be written to the output file breakthrough out followed by the location in X Y and Z space of the point of interest For a reaction path calculation jxplot jyplot and jzplot should be all equal to 1 7 9 Transport Parameters In the last section of the input file threedin the code reads in Cards 20 through 22 information about the transport parameters to be used in the problem The velocities specified are Darcy velocities in units of m m yr The user has the option of specifying either a constant veloc
87. per m mineral Ctot gives the intrinsic surface area of the mineral to be used in computing an initial total number of surface hydroxy sites per mass solution units of m 7 activity of species is fixed by guess j i 8 concentration of species is fixed by guess j i GUESS 1 K initial guess for species concentration In the case of options 7 or 8 guess actually fixes either the concentration or the activity of the species In the case of option 5 primary species is a surface complex guess is the moles surface hydroxyls per m mineral CTOT LK total concentration of species if itype 1 If itype 4 gas constraint ctot is used to specify the partial pressure of the gas g Oz g while guess is used as a guess of the concentration of the basis species related to the gas e g COx2 ag If option 5 primary species is a surface complex is specified then ctot is the intrinsic surface area of the mineral in units of m g7 which will be used in computing a total initial concentration of the surface hydroxyl NCON K name of mineral if itype 3 or gas if itype 4 enclosed in single quotes to be used as equilibrium constraint for named species If itype 5 then this is the mineral upon which the surface hydroxyl sites are developed This ensures that the total number of adsorption sites which otherwise can be treated similarly to other total concentrations in the system is linked to the potentia
88. pre character strings dumlabe1 and exponential de pendences depend DUMLABEL MM character string containing species name affecting far from equilib rium dissolution rate DEPEND MM LL K exponential dependence of the far from equilibrium rate on species given in preceding dumlabel RATEO LL K far from equilibrium mineral dissolution rate constant at 25 C in units of moles m s Can be used as the rate constant at temperature if activation energy is set 0 EA LL K activation energy in units of kcal mole end of mineral read OS3D GIMRT Manual 33 Steefel amp Yabusaki The user needs to be careful in inputting data in this section because of the range of options After read ing the mineral name and nreact the number of parallel reactions involving the mineral the code will read nreact following lines each one corresponding to the parallel reaction Each successive line cor responding to a separate reaction mechanism includes the number of species affecting the far from equi librium dissolution reaction see Equation 28 The code will then read on the same line npre character strings and exponential dependences This is then followed on the same line by a read of the rate at 25 C in units of moles m and the activation energy In this fashion the user can specify multiple parallel reactions e g pH dependent and pH independent or ligand promoted and proton promoted re action mechanisms a
89. ration of the basis species related to the gas e g COa aqy If option 5 primary species is a surface complex is specified then ctot is the intrinsic surface area of the mineral in units of m which will be used in computing a total initial concentration of the surface hydroxy NCON K name of mineral if it ype 3 or gas if itype 4 enclosed in single quotes to be used as equilibrium constraint for named species If itype 5 then this is the mineral upon which the surface hydroxyl sites are developed This ensures that the to tal number of adsorption sites which otherwise can be treated similarly to other total concentrations in the system is linked to the potentially changing surface area of the mineral OS3D GIMRT Manual 55 Steefel amp Yabusaki e NAMK VOLIN LK AREAIN LK NAMK D character string containing mineral name VOLIN LK initial volume fraction of mineral I in the Kth zone AREAIN I K initial reactive surface area of mineral I in the Kth zone Units of Me td Meek e Card 16 NDISTO IXXLO JXXHI D JY YLO JY YHI JZZLO JZZHMD NDIST I geochemical condition number according to the order given in Card 15 to be distributed in space If the geochemical number is negative see second example above the code will treat this domain as fixed throughout the simulation If the geochemical condition is posi tive then the code will initialize the domain to the appropriate geo
90. reas and porosity if desired are up dated The individual mineral volume fractions dm are updated with the expression m t At m t VnTm t At At 73 where is molar volume of the mineral and rm is its reaction rate The user may either choose to leave the porosity constant in which case its value is determined by the user in the input file or the user may choose to update the porosity Where updated the porosity is computed from the 1 S dm 74 Note that in this formulation both mineral volume fractions and reactive surfaces are assumed constant for any one time step This decoupled two step approach normally works well since the volume frac tions of the minerals usually evolve much more slowly than the solute concentrations However in cer tain cases the approach may lead to inaccurate solutions The code has a check to make sure that the volume fractions of a dissolving phase does not go below zero in which case a mass balance error would occur Where the code indicates that this is occurring frequently the time step should be reduced In ad dition one can specify a maximum rate of increase using the parameter vpptmax described below or decrease using the parameter vdissmax described below of the mineral volume fractions The two step approach employed here may require smaller time steps when either the dissolving or precipitating phase has a high solubility or where very little of it is present Where only
91. rential equation for the total concentration of each of the primary species In a 2D problem one may have dependencies of course on other neighboring grid cells in the Y direc tion i e north and south of the grid cell of interest In the case of the operator splitting or sequential iteration approach employed in OS3D the Newton step only involves homogeneous and heterogeneous reactions within the grid cell itself since the transport is handled separately i e there are no nonlinear terms involving spatial derivatives The Jacobian matrix then takes the form se 967 751 C gt 60 The individual terms f Cy make up the elements of the Jacobian submatrices From Equation 55 it follows that the Jacobian submatrices can be written further as Of P A IRT an P P P acy 96 At C 7 OS3D GIMRT Manual 5 Steefel amp Yabusaki ar 80 gt ace om and 9 guy fy y a lt 63 ach ack From Equation 12 and Equation 9 the partial derivatives of the total concentrations with respect to the primary species concentrations are Lichtner 1992 Nz Var dl din y U i Uy Ca e hs be o T aol e 6 06 T Nz din y diny 64 J 1 UC a E 807 Din Ta TE Fj where J is the ionic strength the subscript 2 refers to the secondary species and k is a second dummy subscript referring to the primary species The partial derivativ
92. rid point in Z direction specifying lower left hand corner of a rectangular do main to be initialized with geochemical condition ndist JZZHI 1 grid point in Z direction specifying upper right hand corner of a rectangular do main to be initialized with geochemical condition ndist Card 17 JBD geochemical condition to be used at the boundaries of the 3D rectangular volume The code expects to read 6 numbers here one for each edge of the cube even for 2D 1D and OD problems A negative geochemical condition number is used only by GIMRT trans forming what would otherwise be a no flux boundary into a Dirichlet boundary see section on Initial and Boundary Conditions ee pages ng Lat TOE AT A EIT EA dela y OS3D GIMRT Manual 39 Steefel amp Yabusaki 7 8 Reaction Path Switch and Plotting Information Card 18 specifies whether the code should be run in reaction path mode or batch mode if no discretization is being used i e nx ny and nz 1 In reaction path mode the mineral volume fractions are not updated since we assume that all reaction products are left behind Card 19 contains information about plotting species concentrations at specific points in space In the example given below Card 19 tells the code to write to the datafile breakthrough out the species H and H 60 at grid point 1 1 1 every 2 timesteps 18 CHOOSE BATCH OR REACTION PATH ONLY APPLIED WHERE NXYZ 1 ipath 1 for reaction path i
93. rid points in the system The rate laws used for mineral precipitation and dissolution in IDREACT are based loosely on transition state theory e g Lasaga 1981 Aagaard and Helgeson 1982 Lasaga 1984 The formulation gives the dependence of the rate on the saturation state of the solution with respect to a particular mineral as a function of the ion activity product Qm defined by No Gac 0 27 351 where the a are the activities of the ions making up the mineral and the vm are the stoichiometric coeffi cients The following is the form for the rate law of crystal growth and dissolution of a mineral Lasaga 1984 Steefel and Van Cappellen 1990 Steefel and Lasaga 1994 n Qm lt i Km where k is either the growth or dissolution rate constant a is the activity of a species raised to an em pirically determined power p which affects the far from equilibrium dissolution rate e g H Am is the surface area of the reacting mineral and M and n are two positive numbers which are also nor mally determined by experiment The bars refer to the absolute value of the quantity and the term sgn log Qm Km gives the sign of the expression negative if the fluid is undersaturated and positive if the fluid is supersaturated with respect to the mineral This formulation ensures that the reaction rate rm has the correct sign when n 1 Km is the equilibrium constant for the mineral water reaction written as the destruction of one mole
94. rmulation A of Kirkner and Reeves 1988 Nz Ne L C gt K7 6 i 1 351 N N M Y gt vimsgn log Q Km Amkm 705 7 31 1 0 m 1 j l j 1 Ne 32 where the 1 5 and y s are the activity coefficients for the secondary and primary species respectively The integrated finite difference method is used to convert the differential operator into an algebraic expression and is discussed in the next section zev OS3D GIMRT Manual 9 Steefel amp Yabusaki 3 2 Operator Splitting Sequential Iteration Approach Unlike GIMRT the code OS3D solves the transport terms and the reaction terms separately Splitting of the transport and reaction steps is presently carried out in OS3D using the Strang scheme in which a half transport step is followed by a full reaction step which is in turn followed by another half transport step Strang 1968 Zysset et al 1994 This form of time splitting involves a three step algorithm which requires no iteration and takes the form ea C IE Cp k 1 Ne Nz 33 followed by a reaction step mas cj Pe pprerest a x R 9 1 Ne 34 which is in turn followed by another 1 2 transport step ont Creacted k k ES eeta k 1y No F Nz 35 At 2 In this scheme OR are the concentrations of the individual species at the end of the first 1 2 trans port step U 8009 are the total concentrations at the end
95. roach 12 3 5 Transport Algorithm in Operator Splitting Approach 13 3 5 Advection ple ea e A Be ASA A Sed Ge en ad 13 3 5 2 Diffusion and Dispersion o 14 3 6 Solving the Equations with Newton s Method 14 3 7 Structure of the Jacobian Matrix 16 3 8 Convergence Criteria o 17 3 9 Updating Mineral Properties 17 4 PROGRAM STRUCTURE 18 4 1 Global Implicit y A A A A AAA LD eh An Ge h 18 3 FEATURES OF THE CODE 3000 t aeaa a 6 lap arn BS A EA ee 19 Del ptas a Be Att e dos whe 2 eas aia e aS ees aes Ss 19 5 2 Boundary Conditions 19 5 3 Choice of Primary and Secondary Species 22 5 4 Redox Reactions 20 066 or e ai a ee 23 Steefel amp Yabusaki 1 OS3D GIMRT Manual 10 11 12 5 51 Acid Base Reactions iia o a id o a ja la o e 24 5 6 Basis Switching dd ses E AAA A NR AAA 24 5 7 Activity Coefficient Model lt lt o o o 24 5 8 Initialization Procedure ee 24 5 9 Temperature Gradients 5 9 64 2 Gene Gee OR A PS ee De ES 25 5 10 Updating Porosity and the Transport Coe
96. s the U s The offdiagonal submatrices in contrast contain only contributions from the total soluble concentrations at the neighboring grid points Note that the global matrix A is sparse i e the ratio of zero to non zero entries in the matrix is large when the number of grid points is large e g greater than 100 even though the submatrices reflecting the coupling between the various species due to reactions may themselves be fairly dense Any method used to solve the set of equations described by Equation 59 which arise in the case of the global implicit method should take into account both the sparseness of the Jacobian matrix A In the code GIMRT the equations are solved now using the WATSOLYV package VanderKwaak et al 1995 In the WATSOLV package one may use either the GMRES algorithm or the CGSTAB algorithm developed by the authors Users should consult the manual prepared by these authors VanderKwaak et al 1995 on the various input options which can be specified At the present time all of these input options are in data statements at the beginning of the driver routine gimrt We have found however that for redox problems and probably any similar problem in which there is a large range of concentrations that the convergence criteria need to be set much smaller than indicated by VanderKwaak et al 1995 in their manual We found that the Newton iterations failed to converge when the tolerance was too loose and a a
97. s for the temperature dependence of the parameters are taken from the EQ3 data base 5 8 Initialization Procedure The codes GIMRT and OS3D initialize the solute concentrations at individual grid points over the do main by carrying out a distribution of species using a number of options The initial concentrations of the primary species may be determined by e carrying out a mass balance on a total concentration e g aluminum e carrying out a charge balance on a species e g CIT OS3D GIMRT Manual 25 Steefel amp Yabusaki e a mineral equilibriorn constraint e g H required to be in equilibrium with calcite e fixing the activity or partial pressure of a gas e g COa 4 e fixing the activity e g pH of a species e fixing the concentration of a species and e fixing the number of moles surface hydroxyl sites per mole mineral and the intrinsic surface area of the mineral i e m g This last option is available in OS3D and GIMRT and it allows one to specify the concentration of surface complexes These are specified by giving a site density mole per m mineral and a surface area of min eral per gram mineral The total concentration of surface hydroxyls developed on the mineral converted to mol kg solvent is then computed from this information and the abundance of the mineral specified in the startup file If a mineral is initially is not present but one would like to consider adsorption on this mineral then
98. s is not physically reasonable since it assumes that the diffusive and or dispersive flux through the boundary does not affect the concentration of the exterior node i e it is fixed This boundary con dition would normally be applied where the physical processes exterior to the domain are such that the diffusive dispersive flux through the boundary is not enough to alter the concentrations there One exam ple might be a boundary along a fracture within which the flow is so rapid that the diffusive flux through the fracture wall has a negligible effect e g Steefel and Lichtner 1994 Another example might be a boundary corresponding to an air water interface where a gas like O2 or CO is present in such large abun dance that it effectively buffers the concentrations in the aqueous phase Another example which is often mentioned in this regard is the case where mineral equilibria fixes the concentrations This is dangerous however since the system is only completely determined in the special case where the phase equilibria results in an invariant point i e no thermodynamic degrees of freedom exist This is extremely rare in open systems despite the frequent invocation of this condition in metamorphic petrology Normally since the concentrations will not be uniquely determined a rigorous treatment would require that a speciation calculation with equilibrium constraints be carried out to actually determine what the concentrations at the exterior
99. scale Rate constants at temperature can be computed by providing the appropriate activation energy and the rate at 25 C or by inputting the rate constant at temperature and specifying a zero activation energy OS3D GIMRT Manual 7 Steefel amp Yabusaki 2 4 Example involving Albite Dissolution In this section we give an example involving albite dissolution to show how both the aqueous complexa tion reactions and the mineral reaction rates result in a nonlinear system Here we consider only the ODEs which result from the reaction term i e the equations would be those in the reaction module when using either operator splitting or sequential iteration If a global implicit approach were used one would also have terms involving spatial derivatives in the equation If we consider a geochemical system involving the primary species H20 H AI 3 SiOo qq and Na and the secondary species all in equilibrium with the primary species OH AIOH H3SiOz and NaCl qq then the differential equations take the form 3 d Nat NaCl ANat 2Alt3 Sios 7 de CO Aaw as kaw 1 3 Ten H dt al alb gt TE 3 H 3 d Alt 3 4042 ANat ALO Sion _ l dt Aap a Ka 1 al en Kip gt H i 3 e Nat 0436 AT OH AMLO IO a al hy E Xa 4 a where we have neglected the mass balance of H20 The actual form of the equations solved by both OS3D and GIMRT eliminates the secondar
100. son and Kirkham 1974 Fluid viscosity is calculated from the data presented by Bruges and others 1966 A temperature dependence is also included in the calculation of the reaction rate constants Equation 29 the diffusion coefficients and the fluid densities Reaction rate constants at temperature are obtained nor mally by extrapolating the rate constants for 25 C given in the input file threedin If one would rather input directly the desired rate constant at temperature than the simplest procedure short of changing the source code is to input the rate as if it were the 25 C data while changing the activation energy to 0 kcal mole There is also an option for directly inputting the diffusion coefficient at temperature There are two options for temperature presently implemented in GIMRT The first is simply to run the simulation isothermally either at 25 C or at temperature The second is to specify a linear temperature gradient which operates only in the X direction Many other ways of handling temperature are possible including allowing for a transient temperature field but this must be provided by the user programmer at the present time This can be done with a simple read of temperatures making sure that the dimensions Yay die MT PS me gt Steefel amp Yabusaki 26 OS3D GIMRT Manual of the temperature field match those specified in the code The positions in the code
101. ssolution which can occur within a single time step Useful where either the dissolving phase is very soluble or where very small amounts are present Where very small amounts of a mineral are present inaccurate mineral front propagation may occur where most or all of the mineral present within a grid volume dissolves in a single time step VPPTMAX maximum rate of increase of a mineral s volume fraction units of yr e Card 5 JTEMP TINIT TGRAD JTEMP flag controlling temperature option 0 isothermal 1 linear temperature gradient TINIT temperature in C at grid point 1 TGRAD temperature gradient in C m7 A positive gradient causes the temperature to increase from grid point 1 a negative gradient causes it to decrease e Card 6 MASTER ISWITCH IGAMMA MASTER string up to 18 characters enclosed in single quotes specifying master variable used in controlling time step ISWITCH flag controlling basis switching OO E lt a e ES E lt a k OS3D GIMRT Manual 53 Steefel amp Yabusaki O no basis switching solves for logarithms of concentrations jas not allowed at this time 2 basis switching every Newton iteration Improved method over earlier versions which involves solving for the linear concentrations and using relaxation to prevent overly large Newton steps or negative concentrations Results in reduced number of cases where non convergence occurs especially in redox e
102. t is with the production or destruction of one mole of the secondary species the reactions become E 4 1 ms 081 T Laraga a OH OS3D GIMRT Manual 5 Steefel amp Yabusaki which implies that the total concentrations are Un o CH0 4C aom CoH 20 O aires Chtt Caon 22 One can think of the total concentration then as a linear combination of the concentration of the primary or basis species and the amount of that primary species present in the secondary species The U s there fore account for all of the solute mass in the system even though in certain cases e g at high pH where the concentration of OH is greater than the concentration of H the total concentration of H may be negative As noted above the choice of H2O Ht and Al as the primary species is not unique and it is therefore possible to choose a different set of primary species which is mathematically equivalent If we choose to make Al OH the primary species then the stoichiometric reaction matrix becomes H20 1 A 5 a alee 23 Al OH so that the total concentrations now become Umo Cm0 4C Cop 24 Uy Cy 4C ai t Cop 25 Usong Cayonyz Cait 26 Writing the chemical reactions differently which mathematically is referred to as a change of basis therefore changes the total concentrations of H20 and H even though the total concentration of the Al bearing species is the same In th
103. tant velocity field Code will attempt to read these numbers even if ivelocity 1 but the information will not be used e VXFILE VYFILE VZFILE character strings containing file names with velocity components This field should include 3 character strings even if the code is to be run with a constant velocity otherwise a read error will result e Card 21 ALFL ALFT longitudinal and transverse dispersivities in meters respectively e IDIFFUS DCOEFF PORPOW IDIFFUS flag controlling choice of a diffusion coefficient O default gt 1079 m s7 and an activation energy of 5 kcal mole 1 user provided uses dcoeff DCOEFF user provided diffusion coefficient Do in units of m s only used if idiffus 1 PORPOW the cementation exponent in Archie s Law Equation 5 mr eT ES
104. the concen tration of a basis species especially in the case of redox may drop to a very small value If the logarithms of the concentrations are used it is still possible to solve such problems without basis switching although Steefel amp Yabusaki 16 OS3D GIMRT Manual in redox systems the code sometimes reduces the time steps many times in order to obtain convergence In both OS3D and GIMRT the logarithms of the concentrations are solved for if no basis switching is used and the linear concentrations are used if full basis switching every Newton iteration is used The user should normally begin by running the codes without basis switching since that will be faster if there are no non convergences of the Newton iteration If the code is either not converging within the maximum number of Newton iterations then full basis switching should be tried The user may also want to com pare the iteration count with and without basis switching since the bottom line is the total CPU required to solve a particular problem 3 7 Structure of the Jacobian Matrix It is instructive to examine the structure of the Jacobian matrix which arises from the discretization of the coupled reaction transport equations since it has a bearing on the choice of a method to solve the set of linear equations given by Equation 32 The Jacobian matrices exhibit a block tridiagonal structure in a one dimensional system The blocks making up the global matrix A are N by
105. the total concentrations the U s in any one iteration and then follow this with a distribution of species calculation to determine the concentrations of the individual species see discussion by Kirkner and Reeves 1988 The method implemented here however consists of solving simultaneously for the complexation which are assumed to be at equilibrium the heteroge neous reactions and transport terms This means that the primary species the C s rather than the total concentrations are the unknowns Following the notation of Lichtner 1992 we can introduce the differ ential operator L U Ve u DV U 30 if we drop the density and mass fraction of H20 terms for simplicity and write the governing differential equations in terms of the total component concentrations as LU RF 9 1 No 31 where the reaction term includes only irreversible in our case heterogeneous reactions Note that in this formulation we assume that the diffusion coefficients are the same for all of the aqueous species both primary and secondary This makes the solution of the reaction transport equation much simpler because the secondary species do not have to be individually transported We proceed by substituting Equation 9 into Equation 11 to obtain an expression for the total soluble con centration in terms of the primary species alone Substituting this result into Equation 31 gives essentially a kinetic version of the fo
106. ures that D gt 0 as 0 0 It should also ensure that the porosity will not go negative if the chemical reactions are driven by mass transport since the porous medium locally becomes a closed system 6 INSTALLATION The codes OS3D and GIMRT are designed at the present time to be installed and run on UNIX systems The code has been tested to date on Sun Digital ULTRIX Silicon Graphics and IBM RS 6000 series UNIX based workstations A number of suggestions follow which are designed to make it easier to cus tomize the installation The directories GIMRT and OS3D contain the subdirectories e src where the source files Makefile and include files are located e database where the default data bases master25 data and mastertemp data are lo cated e doc where the documentation in the form of a PostScript file is located The user needs to modify the source code so that the correct address of the GIMRT or OS3D directory tree can be found The address of the data bases is specified in the file database f in the lines open appropriate data base file if datal ne then open iunit5 file datal status ol1d err 334 ed AS mr me lt wee gt 3 2 A 12 A nee e E OS3D GIMRT Manual 27 Steefel amp Yabusaki write gt using database a datal else if temp eq 25 d0 then open iunit5 file home steefel GIMRT database master25 data status o1d err 334 write
107. xamples Does require some extra CPU time however IGAMMA flag controlling method for computing activity coefficients 0 unit activity coefficients except for H20 2 extended Debye Hiickel without calculation of the derivatives of the activity co efficients with respect to the primary species e Card 7 NSTOP PRTINT NSTOP number of graphics files to be printed out 0 no time stepping distribution of species only gt 1 number of output files at prtint intervals PRTINT times at which output should be written at least nstop times should be listed e Card 8 PRIMARY SPECIES e NAMC J name of primary species to be included enclosed in single quotes 1 marks end of primary species read e Card 9 SECONDARY SPECIES e NAMCX D name enclosed in single quotes of the species to be included yi marks end of secondary species read e Card 10 NAMK K NREACT K SATI K SAT2 K THRESH K NAMK K character string containing mineral name NREACT K number of parallel reactions involving the mineral The code will read nreact fol lowing lines SAT1 K order of reaction rate law n in Equation 32 SAT2 K order of reaction rate law M in Equation 32 THRESH K supersaturation log Q K necessary for precipitation e NPRE LL K DUMLABEL MM DEPEND MM LL K RATEO LL K EA LL K NPRE LL K number of species affecting the far from equilibrium dissolution rate The code will t
108. y species entirely writing each in terms of the primary species using Equation 9 3 NUMERICAL FORMULATION GIMRT uses a global implicit or one step method to solve the reaction transport equation In contrast OS3D uses either classic time splitting of the transport and reaction terms or it can be run with the sequen tial iteration method proposed by Yeh and Tripathi 1989 There are some advantages to each approach Yeh and Tripathi 1989 cite the major advantage of the operator splitting or sequential iteration approach as the lower memory requirements of these methods compared to the global implicit and the greater speed with which a single time step can be completed Although these may be advantages of the time splitting approach over the global implicit probably the most significant advantage of the time splitting approach is the ability to use algorithms for high Peclet number transport i e advection dominating over disper sion and molecular diffusion which have less numerical dispersion than those readily usable in a global implicit scheme GIMRT for example using a first order accurate upwind scheme which along with the backwards differentiation leads to a significant amount of numerical dispersion In contrast OS3D uses a third order accurate total variation diminishing or TVD method proposed by Gupta 1991 The TVD algorithm results in signficantly less numerical dispersion than the upwind scheme used in the global implicit i
109. z id Eno In the integrated finite difference formulation a control volume is defined and the differential equations are converted to a set of algebraic equations which include the fluxes across the boundaries of the volume and a source or reaction term Equation 41 can be written Patankar 1980 U Ou e min 1 221 221 f R dz lt 0 42 where the subscripts e and w are the derivatives evaluated at the two boundaries of the control volume Equation 42 can be expressed as a set of algebraic equations using piecewise linear representations of OS3D GIMRT Manual 11 Steefel amp Yabusaki the derivatives Patankar 1980 De u uP Dyluj R Az 0 43 da a where De and D refer to the appropriate value for the diffusion coefficients at the boundary of the volume averaged between the two grid points and the U s superscripts E W and P are the total concentra tions evaluated at the nodal points themselves Figure 1 The term Rj is a function of the primary species at the nodal point P The spacings dx e dx and Az are shown in Figure 1 Ax Control Volume Figure 1 Discretization scheme for a one dimensional problem using inte grated finite differences Following the terminology of Patankar 1980 the discretized equations can be recast in the form apU 05 awU d 0 44 where 45 E dx Dy dz ap QE aw 47 dan

Download Pdf Manuals

Related Search

os3D/GIMRT

os3D/GIMRT

Contents

Download Pdf Manuals

Related Search

Related Contents