Home

a copy of the CrunchFlow manual

Contents

1. V 90 M o C t vX i l N N Velup M o C 9 v X DVIp Myuo C vX i l yf R j 21 N 11 Note that only the term RU remains on the right hand side of Equation 11 because we have assumed that they are the only kinetic reactions Example of a Total Concentration If a total concentration U 3 is defined Reed 1982 Lichtner 1985 Kirkner and Reeves 1988 N U C 2 v X 12 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 O min PLA sM noU V e up M 4U DV p M U R Q L N 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 total concentrations may take on negative values Such quantities however do appear occasionally in geochemistry the best example of which is alkalinity In fact the alkalinity which may take on negative values is just the negative of the total H concentration where CO aq or H CO is chosen as the basis species for the carbonate system Where the system includes sorbed species 1 e surface complexes which are immobile then it is a
2. gt V mbn 7 m 1 where r is the rate of precipitation or dissolution of mineral m per unit volume rock v is the number of moles of 2 in mineral m or stoichiometric coefficient if the reaction is written in terms of the dissolution of one mole of the mineral and N is the number of minerals present in the rock Following convention we take r as positive for precipitation and negative for dissolution Equation 1 provides a general formulation for the conservation of solute mass which makes no assumptions of chemical equilibrium If we assume that the various aqueous species are in chemical equilibrium however it is possible to reduce the number of independent concentrations that is the number that actually need to be solved for Mathematically this means that in a system containing N aqueous species the number of independent chemical components in the system N is reduced from the total number of species by the N 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 N primary or basis 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 N A
3. CH20 label aerobic respiration type monod rate 25C 9 699 activation 0 0 monod terms O2 aq 15 0e 06 inhibition The inhibition terms take a similar if inverted hyperbolic mathematical form K C Ia In i K nm b where Kin is the inhibition constant and Ci refers to the inhibiting species An example involving denitrification 1s given in the following database input CH20 label denitrification type monod rate 25C 10 00 Regnier and Steefel 1999 activation 0 0 kcal mole monod terms NO3 45 0e 06 inhibition 02 aq 30 0e 06 Note that both of these forms are typical of those used for single as opposed to dual Monod formulations In these cases the organic carbon is considered to be in excess of any limiting concentration At this point all Monod type reactions are assumed to be irreversible that is no use is made of the equilibrium constants Irreversible Irreversible reactions are the simplest kind in that no use is made of the equilibrium constants The reactions are assumed to take the form Ea R A k exp as n Akn lt 2 Ta and are therefore similar to the TST form except that a dependence on the saturation state is missing Because of the similarity to the TST reactions the input is also similar for example Iron label tce type irreversible rate 25C 6 92 55 CrunchFlow Manual Carl I Steefel activation 1
4. Hindmarsh Keyword followed by true or false which selects the linear solver for one dimensional problems Syntax hindmarsh logical logical is a standard FORTRAN logical true or false Default true Explanation The keyword parameter hindmarsh is followed by true or false or yes or no and is used to select the linear solver for one dimensional problems Setting hindmarsh as true turns on the direct block tridiagonal solver written by Alan Hindmarsh 1977 Setting hindmarsh to false turns on the PETSc solver Lag_activity Keyword followed by true or false which selects when activity coefficients are updated Syntax lag_activity logical logical is a standard FORTRAN logical true or false Default true Explanation The keyword parameter lag_activity when true causes the code to update activity coefficients only at the beginning of the time step outside of the Newton iteration loop In this scheme activity coefficients are based on the ionic strength calculated from the previous time step Since at this time ionic strength is not included as an independent unknown in CrunchFlow the method yields accurate results in most cases and shows better convergence behavior than does a scheme in which activity coefficients are updated throughout a particular Newton iteration loop The method may not be completely accurate when large contrasts in ionic strength occur in a problem although even this depends on the time step used Later_inp
5. v A 1 N 8 j l where the A and the A are the chemical formulas of the primary and secondary species respectively and v 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 N X K y IIoc 2 1 N 9 jal CrunchFlow Manual Carl I Steefel 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 the destruction of one mole of the secondary species Equation 8 implies that the rate of production of a primary component 7 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 N R vr 10 if 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
6. Ama A Kaw 1 Ki Ht 3 de Er Pra NT a a wg d H OH 4AIOH H Si0 __4 4 gory h wat Ea el dt alb g a 4 a Ht where we have neglected the mass balance of H O The actual form of the equations solved by GIMRT eliminates the secondary species entirely writing each in terms of the primary species using Equation 9 Numerical Approach 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 sequential 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 1 e advection dominating over dispersion 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 numeric
7. If zone is specified then the code expects the X Y and Z coordinates in the form of a beginning JX and ending JX a beginning JY and an ending JY and a beginning JZ and an ending JZ in each case separated by a hyphen If default is specified instead of zone the code will initialize the entire spatial domain 1 e all grid cells to the value More than one specification of tortuosity can be and normally is provided for example tortuosity 0 1 default tortuosity 0 01 zon 1 20 1 1 tortuosity 0 02 zon 1 1 1 20 1 Zones specified later in the sequence overwrite zones specified earlier in the example above the default specification sets the tortuosity value to 0 1 initially in all of the grid cells while the next tortuosity keyword changes the value to 0 01 in grid cells JX 1 to JX 20 The convention adopted here is that all three dimensions must be specified even in a one dimensional problem The tortuosity is considered here to be a number less than 1 some formulations divide by the tortuosity so that it is always greater than 1 The tortuosity 1 as used in CrunchFlow therefore can be viewed as a parameter which allows an effective diffusion coefficient in porous media Dpy to be calculated from the diffusion coefficient in pure water D and the porosity according to Dy PIX Jy JAT ix Jy JZ D Anisotropy_ratioY Keyword followed by a value defining the anisotropy ratio in the Y d
8. In none is provided a single column format consisting of a single porosity per line with NX varying first then NY and then NZ is assumed PEST CreatePestInstructionFile Keyword followed by true or false which selects whether or not to create a PEST instruction PestExchange ins file Syntax CreatePestInstructionFile ogical Default False Explanation This logical determines whether a PestExchange ins file is created with instructions for PEST to read output on cation exchange from CrunchFlow The format specified within the instruction file will match the format used by CrunchFlow with the cation exchange species listed in the order given the exchange keyword within the PEST keyword block Once created this option should be turned off so as to avoid writing over the PestExchange ins file 87 CrunchFlow Manual Carl I Steefel CreatePestExchangeFile Keyword followed by a file name Syntax CreatePestInstructionFile filename Default PestExchange out Explanation This keyword is used to set the filename to be used for exchange output used in PEST runs This option is useful where multiple input files are being run together thus requiring different filenames for the output Exchange Keyword followed by concentration units followed by a list of exchange species to be written out to a file PestExchange out Syntax Exchange units exchange species list Default None Explanation This keyword indicates the list of
9. 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 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 advective 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 1 e there will 26 CrunchFlow Manual Carl I Steefel be a non zero diffusive and or dispersive flux through the boundary In OS3D however the code always assumes the dispersive 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
10. e g aluminum 28 CrunchFlow Manual Carl I Steefel e carrying out a charge balance on a species e g Cl e amineral equilibrium constraint e g H required to be in equilibrium with calcite e fixing the partial pressure of a gas e g CO 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 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 mineral 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 the total concentration of surface hydroxyls on this mineral is set to a small number 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 below Temperature Gradients Since GIMRT is based on finite difference methods it has no problem handling temperatures which are non constant either in time or in sp
11. species and minerals in an actual reactive transport run need to copy the list of secondary species and or minerals over to the keyword blocks for these categories within the input file and then deselect the database sweep option In this mode the code runs essentially like EQ3 does in terms of loading all possible relevant species given a choice of primary component species This option supersedes whatever is set with the speciate only keyword described below Debye Huckel Keyword followed by true or false which selects model for calculation of activity coefficients Syntax debye huckel logical logical is a standard FORTRAN logical true or false Default true Explanation The keyword parameter debye huckel is followed by true or false or yes or no and is used to select the model to compute activity coefficients Currently the only two options are the extended Debye Huckel model Helgeson or activity coefficients equal to 1 no activity corrections Density module Keyword followed by an option either temperature sodium nitrate sodium chloride or potassium nitrate 35 CrunchFlow Manual Carl I Steefel Syntax density module option Default Temperature Explanation The default is to calculate the fluid density only as a function of temperature Limited options are available to calculate the fluid density as a function of composition with sodium nitrate and sodium chloride solution being the two most common Disso
12. 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 90 CrunchFlow Manual Carl I Steefel diffusion Int Journ Numerical Methods in Engineering 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 513 528 Steefel C I 1992 Coupled fluid flow and chemical 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 coupled 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 680 684 Steefel C I and Lasaga A C 1994 A coupled model for transport of multiple chemical species and kinetic precipitation dissolution reactions with application to reactive flow in single phase hydrothermal systems American Journal of Science 529 592 Steefel C I and Lichtner
13. In this format the number of pairs of stoichiometric coefficients and species names is determined by the preceding value in Number of species in reaction if there is a mismatch a read error will result The length of the Log K array is given by the number following temperature points in the first line of the database file An example is given by the entries for quartz and calcite Quartz 22 6880 1 1 00 SiO2 aq 4 6319 3 9993 3 4734 3 0782 2 7191 2 4378 2 2057 2 0168 60 0843 Calcite 36 9340 3 1 00 H 1 00 Ca 1 00 HCO3 2 2257 1 8487 1 3330 0 7743 0 0999 0 5838 1 3262 2 2154 100 0872 Kinetic Database Entries The kinetic database entries include information on rate formulations and coefficients to be used in a simulation Some of these values can be overwritten in the input file at run time With the exception of the case when generic_rates is specified with a value in the RUNTIME block however there must be an entry in the database corresponding to the particular reaction even if some of the values are not used The beginning of the section on solid liquid phase kinetics in the database is marked by the string Begin mineral kinetics which is not enclosed in quotes Each entry is delimited by a line beginning with a which speeds up the database searching This is followed by a string giving the name of the solid phase Names given in the input file described below must matc
14. Ioa e soouss cider eroista iss AOD Spatial Profile NR TED METTE Time series R Time series print CU Time series interval eseese 46 Time series output eseese 46 IVEAICO NTO VIE AEE eerte b ERE ERE AITA A E E ET Li BILDER IRI ccc ERA i BALDI 2 RERUM ME ERE M PET OMEOS EEREN E EE vi esee eaten Peso u iuo dates vint eau voee E NU dues ete E esos voee tese do e vod eee Caesar ad Database Rormiat 2oieduee o dinero ul stabs edoceri oit Oe oae o Uic od uR oto iss dope wa sevo deasdeas kae deas oi dew 40 Input File Entry of Primary Species 4 eere ecce esee ee ee ee eere eese eee esse sess sese es sees S eeeesess 49 ies Database Ermats oes oce eti tope ore itas e aves vea e i eo uv oai ei L Input File Entry of Secondary Species eene ee ee eee e ee eee eee ee ee eese eese eeesssee OL CrunchFlow Manual Carl I Steefel Database Bornidts 5 sc cusscensssescoseanaconspetccsinivdcciisestisccascvacasevebaacetencsdosieouctosnseassongansesesisiod L Thermodynamic Database Entries aee p rhe bd e se quad aed ae ada 2d Kinetic Database Entries e Urt teet rr Ue Hx URP EET ae 22 Entries in Input File Sees eiecti OO AQUEOUS KINETICS paio ER EE E DI Database HOrmiatsiesiessesescasossesverosnsevesonguus cscudsscesespoonss sonncvselengsotescospanssensbovsvenssevasenssebdsccsecs OU Types of RAL LAWS qc s
15. P Ww D u u D u u dx dx R Ax 0 37 where D 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 concentrations evaluated at the nodal points themselves Figure 1 The term ke is a function of the primary species at the nodal point P Following the terminology of Patankar 1980 the discretized equations can be recast in the form apU a ay d 0 38 where D 39 amc D 40 ay dx 40 Ap a y 41 16 CrunchFlow Manual Carl I Steefel d RM Ax 42 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 38 One way to represent the first derivative which appears in the advection term is to use a central difference method Z Ug Up 43 OX Jp 2Ax where the subscript P indicates that the derivative is evaluated at the nodal point P and Ax is 1 2 6x x The truncation error for the central difference representation is said to be of the order Ax because in the Taylor series expansion upon which it is based only quadratic terms and higher are neglected Lapidus and Pinder
16. Species Activity Constraint Specifying a species activity 1s also straightforward and takes the form lt PrimarySpeciesName gt Concentration activity A specification of the pH is a special case which does not require the trailing activity keyword pH 7 8 Gas Constraint The gas constraint is used to calculate a primary species activity so as to obtain equilibrium with respect to a gas at a specified partial pressure fugacity At this point no fugacity corrections are available in the code so effectively partial pressure fugacity Partial pressure is input in bars The format for using the gas constraint is to follow the primary species name by the name of the gas with which the species is to be equilibrated and then the partial pressure in bars primary species name gt gas name gt lt partial pressure bars To adjust the O2 aq activity to be at equilibrium with respect to atmospheric oxygen one would use the following form 02 aq 02 g O20 Mineral Constraint The mineral constraint operates in a similar fashion to the gas constraint except that no trailing value is needed This is because minerals or solid phases are assumed at this stage to be pure that is their activities 1 Unlike some other codes like PHREEQC CrunchFlow does not carry out a formal reaction path calculation to attain equilibrium with respect to the mineral The mineral constraint is used together with the other constraints in the geochemical system
17. Steefel with differing labels and rate laws to capture both the pH independent and pH dependent rates After discussing the remainder of the input for the solid liquid phase kinetics an example of multiple parallel reactions involving calcite will be given Several types of reactions may be specified The five options currently available are type tst type monod type irreversible type PrecipitationOnly type DissolutionOnly and are described briefly below TST These three different reaction types involve different kinds of input following the rate 25C because of the difference in the rate formulations The tst rate law is described in detail by Lasaga 1981 1984 and by Aagaard and Helgeson 1981 and takes the form R A k ep z m a exp m g jr AG g In Q RT K where A is the mineral surface area k is the intrinsic rate constant in units of mol m s E is the activation energy kcal mole Q is the ion activity product for the mineral water reaction Keq is the corresponding equilibrium constant and Le is a product representing the inhibition or catalysis of the reaction by various ions in solution raised to the power n Note that this last term operates far from equilibrium that is when Qm K eq is very small i e the mineral is very undersaturated Following the specification of the activation energy is the specification of the dependences of the far from equilibrium rate on vari
18. and Z components In a one dimensional problem in X the Y and Z velocities are optional Units are changed inside the FLOW block with the space units and time units keywords Pump Keyword followed by a value s giving pumping rate the name of the geochemical condition to be used in the case of an injection well and its grid coordinates Syntax pump rate label jx jy lt jz gt Default None Explanation This keyword beginning with pump sets a pumping rate positive for injection negative for withdrawal in units of l sec and assigns a geochemical condition to it This 1s followed by the coordinates of the well with JY and JZ being optional in the case of a one dimensional problem The geochemical condition is identified by its label although the condition is only used in the case of an injection well Gaspump Keyword followed by a value s giving gas pumping rate the name of the geochemical condition to be used in the case of an injection well and its grid coordinates Syntax gaspump rate label jx jy lt jz gt Default None Explanation This keyword beginning with gaspump sets a gas pumping rate positive for injection negative for withdrawal in units of l sec and assigns a geochemical condition to it This is followed by the coordinates of the well with JY and JZ being optional in the case of a one dimensional problem The geochemical condition is identified by its label although the condition is only used in the c
19. condition Syntax set_porosity value Default Provided by value in POROSITY keyword block Explanation This keyword overrides the value set in the POROSITY keyword block If absent that value is used for the geochemical condition NOTE This option only works if fix_porosity is set in the POROSITY keyword block Otherwise porosity is calculated from the 1 minus the sum of the volume fractions of the minerals Set_saturation Keyword used to set a liquid saturation within an individual geochemical condition Syntax set_saturation value Default Provided by value in RUNTIME keyword block otherwise 1 0 Explanation This keyword overrides the value set in the RUNTIME keyword block If absent that value is used for the geochemical condition Types of Constraints Aqueous Species Various types of constraints may be specified for aqueous species in the distribution of species calculations carried out on each of the geochemical conditions These constraints are summarized in Table 4 Constraint Type Use Example Total Concentration Mole balance on total aqueous or Na 0 001 total aqueous sorbed concentrations Species Specify an individual primary Na 0 001 species Concentration species concentration Species Activity Specify an individual primary Nat 0 001 activity 65 CrunchFlow Manual Carl I Steefel species activity pH Specify activity of hydrogen ion pH 7 8 Gas Equilibrate primary spec
20. every grid cell in the domain is specified implying that this block is read after the DISCRETIZATION keyword block Example NITIAL CONDITIONS quartz zone 1 42 1 42 portlandite zone 227533 1 42 nonreactive zon 12 42 12 42 fix END 71 CrunchFlow Manual Carl I Steefel BOUNDARY_CONDITIONS Boundary conditions are required for any problem involving transport The input currently involves specifying boundary conditions corresponding essentially to ghost cells outside of the domain at the Syntax x begin condition name boundary condition type x end condition name boundary condition type y begin condition name boundary condition type y end condition name boundary condition type z begin condition name boundary condition type z end condition name boundary condition type Default None Explanation Boundary conditions are set for ghost cells outside the faces of the 3D domain Currently uniform boundary conditions are specified over the entire face although this is being revised to provide more flexibility For additional flexibility the user may fix initial conditions at the first or last grid cell in a coordinate direction thus making it possible to specify multiple geochemical conditions on any one face For fixed flux problems it may be easier to make an injection well with the correct flux at the grid cell along the boundary Boundary condition types are Dirichlet use
21. file will be every time step Since this is a 1D problem the Y and Z coordinates need not be given Spatial profile Option to write a spatial profile or snapshot at the specified time s Syntax spatial profile timel time2 time3 amp timen timel time2 etc are real numbers giving the time at which the spatial profile is to be written The ampersand can be used to continue listing of output times on multiple lines Any one line can be up to 132 characters long Default If not specified the code assumes that no output file is to be written and therefore no time stepping will occur Explanation This keyword parameter instructs the code to output a spatial profile or snapshot at the specified time The time units are taken from the time units value specified in the OUTPUT keyword block if absent time units of years are assumed An ampersand can be used to continue listing of output times on the following line This 44 CrunchFlow Manual Carl I Steefel important parameter in addition to specifying when spatial profiles will be written controls whether time stepping is carried out at all The code will run until the last output time in the list is reached unless the parameter steady_state in the RUNTIME keyword block is turned on in this case the code runs until a steady state is achieved or until the final spatial profile time is reached whichever is first Field variables written out as part of the spatial prof
22. files will be created each time a set of spatial profiles are written see keyword spatial profile in the OUTPUT keyword block With a normal successful termination of the code the restart file will be updated at the end of the run Subsequent output to the breakthrough files see keyword time series can be appended by adding the optional append keyword after the filename In the same way the file counter n for example pHn dat will be updated as well if this append option is chose so spatial profiles written after restart will reflect the files written previously before restart Save restart Option to save a restart file Syntax save restart filename filename is an optional string up to 132 characters indicating the name of the restart file to be written to Default crunch rst Explanation This keyword parameter specifies that the code should create a restart file with the name provided If the keyword save restart is not provided restart information will be saved to the default restart file name of crunch rst Screen output Keyword followed by an integer value giving the interval at which run time output is written to the screen Syntax screen output integer value Default Explanation This gives the interval at which run time output time time step size number of Newton iterations required is written to the screen This is useful where the code is taking a large number of time steps to complete a ru
23. 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 the 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 Choice of Primary and Secondary Species After setting up the problem physically one needs to decide on the chemistry to be included GIMRT requires that one specify an initial choice of primary species which then determine the number of independent chemical components in the system This choice of primary species is not unique however as discussed in the section 2 2 The user must choose however 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 ar
24. it involves considerably more CPU time At the present time therefore we 15 CrunchFlow Manual Carl I Steefel recommend using the Strang operator splitting scheme rather than sequential iteration We have incorporated the sequential iteration approach nontheless and describe it below Integrated Finite Difference Formulation The set of partial differential equations represented by Equation 34 is discretized using the integrated finite difference method Patankar 1980 Marsily 1986 As an example of how the discretization proceeds 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 2 OU min a 9D R 0 E 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 35 can be written Patankar 1980 U U Ean pp pese f R dx 0 36 ax j OK j w where the subscripts e and w are the derivatives evaluated at the two boundaries of the control volume Equation 36 can be expressed as a set of algebraic equations using piecewise linear representations of the derivatives Patankar 1980 E P
25. name is the name of the exchanger given in the database Multiple exchangers can be listed If the exchanger name is not followed by anything it is assumed that exchange is on the bulk material in which case a CEC is calculated from the combination of the solid phase density porosity and liquid saturation see Section 1 2 14 Alternatively it is possible to specify exchange on a specific mineral in which case the CEC is calculated from the volume fraction of the mineral which may change with time for example exchange Xkaol on Kaolinite Also required is the keyword convention convention activity convention where activity convention must be either Gaines Thomas of which Gapon is a variant or Vanselow SURFACE COMPLEXATION Surface complexation follows the approach outlined in Dzombak and Morel 1990 with either a double layer or non electrostatic model possible Currently complexation must be on a specific mineral so a valid mineral name listed in the MINERALS keyword block must be given for example gt FeOH strong on Fe 0OH 3 gt FeOH weak on Fe OH 3 To specify a non electrostatic model the mineral name should be followed by the keyword no edl for example 63 CrunchFlow Manual Carl I Steefel gt FeOH strong on Fe OH 3 no edl gt FeOH weak on Fe OH 3 no edl The capability for surface complexation on the bulk material will be added soon GEOCHEMICAL CONDITIONS Geochemical conditions are keyword blo
26. of a secondary phase not initially present one can specify a threshold mineral volume fraction which is used to calculate the bulk surface area until the computed time evolving volume fraction exceeds the threshold value A threshold value for the case of an initial volume fraction 0 is set by name of solid phase 0 0 specific surface area value threshold volume fraction 69 CrunchFlow Manual Carl I Steefel So for example we could specify that the bulk surface area of calcite which is initially not present be calculated from a threshold volume fraction of 0 0001 and a specific surface area of 2 m g with the following form Calcite 0 0 specific surface area 2 0 0 0001 With this formulation the bulk surface area of calcite would be calculated directly from the calcite volume fraction and specific surface area once it exceeded the threshold value of 0 0001 set her This is a simple if not particularly mechanistic way to capture a nucleation event Surface Hydroxyl Concentrations If surface complexation is included in a particular problem then surface hydroxyl concentrations must be specified in each of the geochemical conditions there is currently no default value available The information which needs to be provided is the surface hydroxyl site density Aes in units of moles sites per m mineral This information is combined with the specific surface area A to calculate a concentration of surface hydroxyl site
27. present in the problem In a one dimensional problem however 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 where Ax refers to the grid spacing at any particular point in space Only for Pe To take advantage of the greater accuracy of the central difference method at small values of Pe a power law scheme proposed by Patankar 1980 is used This method slides between a grid 17 CrunchFlow Manual Carl I Steefel fully centered form at Pe lt 2 to an upstream weighted formulation at Pe gt 10 where the upwind scheme is necessary for stability Using the power law scheme of Patankar 1980 the finite difference coefficients become 5 a D 0 1 228 0 v 46 D TS D ay D of 1 248 0 v 47 dp a dy 48 where the operator A B means that the greater of the values A or B is used equivalent to the FORTRAN operator AMAXI A B and A refers to the absolute magnitude of the value A fuller description 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 le
28. selected with the gravity keyword and appropriate boundary conditions typically just a fixed pressure at one end of the domain TEMPERATURE This keyword block is used to set various temperature parameters A full calculation of the temperature will be added soon but in the meantime temperature can be set block by block using individual geochemical conditions or by reading temperature values from a file Set_temperature Keyword followed by a value giving the default temperature Syntax set_temperature value Default 25 C Explanation The keyword parameter set_temperature is used to set a global value for the temperature The value provided here can be overwritten by specifications in individual geochemical conditions provided these conditions are distributed in space as initial conditions Temperature_gradient Keyword followed by a value giving the temperature gradient Syntax temperature gradient value Default 0 0 Explanation The keyword parameter temperature gradient can be used to set a temperature gradient in units of C m A positive value means that the temperature gradient will increase from grid point 1 to grid point NX The temperature gradient when input in this way operates only in the X direction Read TemperatureFile Option to specify a file with temperatures to be read Syntax read temperaturefile filename format 85 CrunchFlow Manual Carl I Steefel filename gives the name of the file up to 132 cha
29. species participating in exchange reactions e g Nat Ca Mg to be written out to the file PestExchange out which can then be used as data or observations in a PEST run using CrunchFlow If the keyword CreatePestInstructionFile is also set as true then the format for the exchange species output will be written to a PEST instruction file PestExchange ins The exchange species list is preceded by a specification of the units which may include Units Alternate Designation mol g mole g moles g mmol g mmole g mmoles g umol g umole g umoles g equiv g eq g equivalents g mequiv g milliequiv g meq g milliequivalents g uequiv g microequiv g ueq g microequivalents g EROSION Read_BurialFile Keyword followed by a file name containing erosion burial values over the X direction Syntax read_BurialFile filename format Backwards compatibility read_burial Default None Explanation This keyword provides a name for a file containing burial erosion rates values in the X coordinate direction Other coordinates Y and Z are not currently supported Units are determined by space_units and time_units specifications within the EROSION keyword block Default units are m yr The format of the file depends on the format specified after the file name but the values of the erosion burial rates are assumed to be listed with NX varying first then NY and then 88 CrunchFlow Manual Carl I Ste
30. terms takes the form lt maximum rate gt lt number of Monod terms integer gt lt electron acceptor or donor gt lt half saturation constant gt lt electron acceptor or donor lt half saturation constant gt where the electron donor or acceptor is a species name which again can be prepended with the string tot to designate a total concentration and the half saturation constant is the concentration of the electron acceptor or donor at which the rate is 2 of its maximum value Inhibition terms of the following form may also be included K C I In K nm where Kin is the inhibition constant and C is the concentration of the electron acceptor or donor Inhibition terms in the aqueous kinetics section are specified by the string Inhibition enclosed in single quotes followed on the same line by Number of inhibiting species lt SpeciesName gt lt Inhibition constant lt SpeciesName gt Inhibition constant Irreversible Irreversible reactions are the simplest kind in that no use is made of the equilibrium constants The reactions are assumed to take the form R kJa and are therefore similar to the TST form except that a dependence on the saturation state is missing Such a reaction rate law can be used to an exponential decay reaction In combination 61 CrunchFlow Manual Carl I Steefel with a series of other decay reactions it can be used to describe an entire decay chain The de
31. the following form Appelo and Postma 1993 written here assuming the Cs is the relevant cation of interest Cs V MX i o CsX i ym where M is the competing cation Nat K Ca mis its charge and X i refers to the ith type of exchange site In the Gaines Thomas convention each exchange site X i has a charge of 1 The activities of adsorbed species correspond to the charge equivalent fractions i fa cr n XO 62 CrunchFlow Manual Carl I Steefel where zy is the charge of cation M q i m is the concentration of adsorbed cation M in exchange site i moles g and the square brackets denote activities The Gapon activity convention is obtained by writing the reactions in every case with a single exchanger Appelo and Postma 1993 Alternatively the Vanselow convention Vanselow 1932 describes the exchanger activity with mole fractions BOy XO The exchange reactions can then be used to write a mass action equation for binary Cs M exchange BO ICs _ XO ICs M Cs BG M X DM T In a single site ion exchange model the CEC is equal to the sum of the charge equivalent concentrations of the adsorbed cations CEC gt u Mau while in a multi site model the CEC is the charge summed over all of the cation exchange sites Cernik et al 1996 Voegelin et al 2000 CEC 5 Sv In CrunchFlow exchangers may be specified with the format exchange exchanger name where exchanger
32. to provide the same number of constraint equations as there are unknown and this system of equations is solved without any path dependence The format for a mineral constraint is simply lt primary species name gt lt mineral name gt Charge Constraint The charge constraint may be somewhat more difficult to implement because the user may choose to balance on an anion when the remainder of the solution is already negatively charged or vice versa In this case without the proper safeguards the concentration of the charge 67 CrunchFlow Manual Carl I Steefel balancing species would go to 0 If this situation arises an error message will result and the run will terminate The format for specifying a charge balance is primary species name charge Solid Phase Inputs In most cases information on mineral concentrations or volume fractions and mineral surface areas are needed to completely describe a geochemical condition Those cases which don t require such information are typically boundary conditions or source terms where only the aqueous phase is entering the domain In this case it is possible to leave the mineral input out of the condition The code does not check to see whether mineral input is or is not needed f it doesn t find any information on a mineral or solid phase which has been loaded in the MINERALS keyword block it will assume its volume fraction 0 and use a default bulk surface area of 100 m mineral m por
33. to the macroscopic transport rates The rate laws used for mineral precipitation and dissolution are based loosely on transition state theory e g Lasaga 1981 Aagaard and Helgeson 1982 Lasaga 1984 This 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 Q defined by N Q a D j 1 where the a are the activities of the primary species used in writing the dissolution reaction for the mineral In order to incorporate the strong pH dependence of most mineral dissolution and precipitation reactions far from equilibrium parallel rate laws are used which are summed to give the overall reaction rate law for a particular mineral Ns dk l 1 ECAN 1 2 28 where k is the far from equilibrium dissolution rate constant for the th parallel reaction p N N Pi IE i l r A gives the exponential dependence on species i of the th parallel reaction 1 e the reaction 11 CrunchFlow Manual Carl I Steefel order K is the equilibrium constant N is the number of parallel reactions and A refers to the surface area of individual minerals in the rock matrix or fracture The exponents n and M allow for nonlinear dependencies on the affinity term and are normally taken from experimental N N studies The term IL a incorporates the effects of various ions in sol
34. with Newton s method which 19 CrunchFlow Manual Carl I Steefel makes use of a Taylor series expansion to linearize the equations Using Newton s method the linearized set of equations for the 1D global implicit case become N of N of N of 7 2 acre 2 acres 2 acevo F 60 iN Jj E Jj z J where f is the discretized differential 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 direction 1 e north and south of the grid cell of interest The individual terms Of oC make up the elements of the Jacobian submatrices From Equation 49 it follows that the Jacobian submatrices can be written further as of OU r oR a Ax 51 OC ac at ac GD of wou A 52 20 act o2 and of our 53 oc oC 92 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 of oC as forming the elements of the global Jacobian matrix A then Equation 50 can be written in the familiar form A C b 54 Once the C are computed they are used to update the concentrations of the primary species c cm C 7 1 N xM 55 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 c
35. 1 00 Currently it is only possible to specify a dependence on reaction affinity in the database Only those rate laws in which a dependence on reaction affinity occurs TST PrecipitationOnly and DissolutionOnly will be affected Monod Monod reactions take a slightly different form than do TST type reactions Following the specification of the activation energy the Monod option has a field for inputting various Monod terms that is the dependence of the reaction rate on electron acceptors and or electron donors The specification of Monod terms takes the form monod terms electron donor or acceptor lt half saturation constant gt electron donor or acceptor lt half saturation constant The half saturation constant is assumed to be in concentration units of moles kg water This in turn is followed by a field for the specification of inhibition terms also in moles kgw inhibition terms inhibitor species name gt inhibition constant inhibitor species name gt inhibition constant The Monod terms take the standard form 54 CrunchFlow Manual Carl I Steefel R ae nl Cow K uu with the quantities in parentheses representing the Monod term Multiple Monod term may be specified although the most common approach is to use a dual Monod form in which a dependence on an electron acceptor and on an electron donor are provided An example involving aerobic respiration is given by
36. 1982 In contrast a forward difference representation given by Z n O A 44 Ox Jp Ox 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 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 relative importance of the two can be described by the dimensionless grid Peclet number v AX Pe 45 D SD grid 2 is the central difference 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
37. 5 0 kcal mole dependence TCE aq 1 0 where the rate of dissolution of zero valent iron is here considered to have a first order dependence on the concentration of TCE in solution PrecipitationOnly Invoking this option in the database causes the code to calculate mineral precipitation only dissolution is suppressed An example is given by KaoliniteYang label Yang type PrecipitationOnly rate 25C 13 47 activation 0 0 kcal mole dependence AffinityDependence 1 00 2 07468 1 00000 DissolutionOnly Invoking this option in the database causes the code to calculate mineral dissolution only precipitation is suppressed An example is given by KaoliniteYang label Yang type DissolutionOnly rate 25C 12 9393 activation 0 0 kcal mole dependence AffinityDependence 1 00 0 5 1 00 This option is useful for specifying different rate laws for precipitation and dissolution where both PrecipitationOnly and DissolutionOnly are invoked separately for the same mineral Entries in Input File The simplest possible form for specification of a solid aqueous phase reaction in the database is to give the solid or mineral name on a separate line within the MINERALS keyword block without any accompanying information If a database sweep has been specified previously then the list of minerals may be copied from the file input file name prefix out to which they have been written and then pasted
38. 6 Zysset A F Stauffer and T Dracos 1994 Modeling of chemically reactive groundwater transport Water Resources Res 2217 2228 9
39. CrunchFlow Software for Modeling Multicomponent Reactive Flow and Transport USER S MANUAL Carl I Steefel Earth Sciences Division Lawrence Berkeley National Laboratory Berkeley CA 94720 USA CISteefel g bl gov October 12 2009 CrunchFlow Manual Carl I Steefel CrunchFlow Manual Carl I Steefel Mathematical Formulatlon oiecsecsaeenserekesnset cce nve sees seisa onpsenssongundsssatssvssespoetes suanevsdengsetescotoucse 2 Conservation Equatiolis 5 p iet US pr NEUEN ERO NM Ue I iN Pe M deste eee iot iain ed Example of a Total Concentration eseeesees eee n ete eere et PER EURN E E iens tapete eee race eur thesis D Reaction Rate Laws eed ET Example Involving Albite Dissolution Diss 12 Numerical Ve cir ree La Global Implicit Approde lviccccisienas pscevosckeroesecnsesuocsdoanssesneseeeasdoavedecusonsovecousesesdesswegsennvesenncnevnees L Operator Splitting Approach oce teta eO Nea Rie reet Esci neve bea c Na ect desine L5 Integrated Finite Difference Formulation e ecce esee ee eee eere ee eee ee ee se seessseeesss LO Transport Algorithm in Global Implicit Approach eeeeeeee eene eren ee eee eees see L7 Transport Algorithm for Operator Splitting Approach s esssisuse LO P tigcorp seis sess TD t eme LS Diffusion and Dispersloll 46e ese ys ee een toa eR Se P ceed n
40. In each case the use of the ppm units refers to the parts per million of the solute using the molecular weight of the primary species to which it refers Therefore if the primary species is NO3 then the concentration in ppm is calculated using the molecular weight of NO3 not for example using N Equilibrate surface When specified on a separate line this keyword instructs the code to partition all of the total concentrations between the aqueous and solid phases i e surface hydroxyl sites and exchangers It does not partition any mass through precipitation dissolution reaction This instruction will only apply to species for which a total concentration constraint is used described below This instruction can also be specified on a species by species basis as 64 CrunchFlow Manual Carl I Steefel discussed below If the equilibrate surface instruction is absent then the code assumes that the total concentration refers only to the aqueous phase and the exchanger and surface hydroxy concentrations are calculated accordingly Temperature Keyword used to set a temperature within an individual geochemical condition Syntax temperature value Default Provided by value in TEMPERATURE keyword block Explanation This keyword overrides the value set in the TEMPERATURE keyword block If absent that value is used for the geochemical condition Set_porosity Keyword used to set a porosity within an individual geochemical
41. MRT false Syntax Courant_number value value is a positive real number between 0 0 and 1 0 Default 0 5 Explanation This is an optional parameter that controls the maximum CFL number by limiting the timestep At when the OS3D option is used The keyword is ignored if the GIMRT option is selected The Courant number cannot be greater than 1 0 for the OS3D option since results in an unstable method 34 CrunchFlow Manual Carl I Steefel Database Option to specify a database file Syntax database filename filename is an optional string up to 132 characters indicating the name of the database file Default datacom dbs Explanation This specification of the name of the database file may also go in the DATABASE keyword block Database_sweep Keyword followed by true or false which selects whether to carry out a full sweep of the database loading all possible species and minerals Syntax database_sweep logical logical is a standard FORTRAN logical true or false Default false Explanation The keyword parameter database_sweep is followed by true or false or yes or no and is used to select when true the option to sweep the entire thermodynamic database so as to load all possible species and minerals In this case the user s choice of secondary species and minerals listed in the input file is ignored After the sweep each of the geochemical conditions is speciated and the code stops User s who want to use these
42. Mineral surface area formulations based on grain size like that presented by Lasaga 1984 cannot be used to describe the mineral surface area in contact with fluid as the porosity goes to zero a porous medium made up of closest packed spheres cannot have a zero porosity Here we use the expression A A 2 G1 o where A and refer to the initial surface area and initial porosity respectively Example Involving Albite Dissolution In this section we give an example involving albite dissolution to show how both the aqueous complexation 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 12 CrunchFlow Manual Carl I Steefel equation If we consider a geochemical system involving the primary species H O H Al SiO 5449 gt AIOH H SiO and NaCl and Na and the secondary species all in equilibrium with the primary species OH aq then the differential equations take the form Fe 3 d Na NaCl e A a k A 1 eT A Sia Je s dt alb H alb a a H d Sio H SiO I actu T 4 3A ar Kay 1 E Je gt Ht 3 d Al AIOH A pat ars Sig an di
43. P C 1994 Diffusion and reaction in rock matrix bordering a hyperalkaline fluid filled fracture Geochim Cosmochim Acta 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 26577 2677 Strang G 1968 On the construction and comparison of difference schemes SIAM J Numer Anal 506 517 Valocchi Albert J and Malmstead Michael 1992 Accuracy of operator splitting for advection dispersion reaction problems Water Resources Res bf 28 1471 VanderKwaak J E Forsyth P A MacQuarrie K T B and Sudicky E A 1995 WATSOLV Sparse Matrix 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 University Press Cambridge 176 p Walter A L Frind E O Blowes D W Ptacek C J and Molson J W 1994 Modeling of multicomponent reactive transport in groundwater 1 Model development and evaluation Water Resources Res 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 11
44. REAR RUNE aose EE Re EXER CO ERROR soser essais SO CrunchFlow Manual Carl I Steefel Fix SAturatlOn ireeie seco eee ceo eed oo eee otc ea va oo cesescssscvsedss susscsssscevetsesdtesssescsvedetestscssetcevstsosstess OO Generi CIS TRECEEEEOUEE EesTr Gimrt pc NIS A Gimrt pclevel rei I Gimrt solver ssi OT Graphics Aeros Hindmarsh esses sentes 22 ec tena eae eoo sesoses ne baee adea bec eene Peres ao ee coc sone b sos e oe Cocos so po Te ue ka c Coc oc so bec oe eade cago UD Lag activity er M S I Later inputhiles oie REPRE E IORSNNERNNYE NO RR EE RH OO Master 1 sees eeee ee ee orae inb adeo eo ee beoe dueta dso cas odotaa eoe or aee No o eese oodo eed eoo eoe ees u oseese rasana ese eo sene OP OVvershoot POleraniee de MEDIIS M2 D mec e ena aU Keaction palthicoi diua aod LC reir Tasa b VEO F HD EE AREE EDe FER LL ENCARTA rEdiu tenem tienen aU Read SaturationFile 40 Restart eere 41 Save restart ut 41 Screen output dias SOLVER ndn inpr qe p Rep pU FEMPCR erpr FOR VR E EQUL QUEUE D NUR De o EAR P Specrite Only 5o ps DR o ERI II ncn Rau Mna dM Steady Stat E MPa Tine Step HANAN cies 200 en hee aditu AE CL oU MO MR Ro D Ce Ta a DER oEE La Te da npE M Rada R Timestep MIELE TERT eH RTT E Time tolerants niic eite eii cett Los poni ae bud dra D beau Cosa a Rd de Qva
45. a e g from a kinetic experiment with output from the simulation MakeMovie Option to write out a movie file of specified total concentrations Syntax MakeMovie speciesl species2 where species and species2 are primary species Default No movie file will be written Explanation This keyword parameter provides the option to write out a total concentration for a primary species in 1D or 2D The interval at which data is written out is given by time series interval DISCRETIZATION This keyword block is used to set up the structured grid for a particular problem The code assumes presently that rectilinear coordinates are used The three keyword parameters which are recognized are xzones yzones zzones each being used to set the grid dimensions in one coordinate direction Xzones Discretization in the first X direction Syntax xzones cells spacing fcells spacing amp cells spacing where Zcells 1s an integer and spacing is a real number 46 CrunchFlow Manual Carl I Steefel Default 1 0 if xzones is not specified Explanation This keyword parameter specifies the discretization in the X first coordinate direction In one dimensional problems the X coordinate rather than the Y or Z coordinate should be used The discretization is given by specifying pairs of integers and real numbers which give the number of cells of a given spacing The following example would result in 10 cells of 1 meter spacing 10 ce
46. ace If either a non default thermodynamic 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 data or mastertemp data In addition the algorithm for fluid density must be modified as well Presently the fluid density is calculated as a function of temperature from a polynomial fit to the data in Helgeson 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 30 the diffusion coefficients and the fluid densities Reaction rate constants at temperature are obtained normally by extrapolating the rate constants for 25 C given in the input file 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
47. al bulk surface area specified is used as long as precipitation occurs if this phase later dissolves the above formulation is used but with an arbritrary initial volume fraction of 0 01 Clearly the expressions above are only a very approximate means of calculating the reduction in surface area of a secondary phase if later dissolution occurs More accurate for secondary minerals and therefore recommended is the use of specific surface area Specific surface area When specific surface area is selected in the input file the reactive surface area used for mineral precipitation and dissolution is calculated from A peci MW Apuk ee y This is the preferred approach for precipitation of secondary phases in particular since the bulk surface area BSA approach currently uses the initial surface area only 1 e the reactive surface area of a precipitating phase does not change as the secondary mineral concentration evolves In contrast with the specific surface area option evolution of the mineral volume fraction causes the bulk surface area to evolve according to the equation above initial If the initial volume fraction of the solid phase is gt 0 then the A77 is calculated directly from the specific surface area If the initial volume fraction 0 that is we are dealing with a secondary mineral initially not present at all the specific surface area cannot be evaluated without some other procedure To handle this case
48. al dispersion In contrast OS3D uses a third order accurate total variation diminishing or TVD method proposed by Gupta et al 1991 The TVD algorithm results in signficantly less numerical dispersion than the upwind scheme used in the global implicit 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 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 13 CrunchFlow Manual Carl I Steefel 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 metho
49. arameters specified by zone e g pressure tortuosity See the keyword permeability_x for details on how this keyword is used to specify permeabilities by zone The convention adopted here is that all three dimensions must be specified even in a one dimensional problem Gravity Keyword used to indicate the coordinate direction of the gravity vector Syntax gravity lt AngletoX gt lt Angle to Y gt lt Angle to Z gt direction Default 90 Explanation This keyword provides the angle relative to the X Y and Z coordinate directions for the gravity vector The influence of gravity is then calculated from g cos where 0 is the angle of the gravity vector relative to the specific coordinate direction X Y or Z Direction can be either down or up and is used to indicate whether the direction the gravity vector is assumed to operate for example toward X NX or toward X 0 So for example the following input gravity 90 0 0 0 90 0 down would be used to set the gravity vector in the Y direction with gravity oriented down toward JY 0 83 CrunchFlow Manual Carl I Steefel Pressure Keyword used to specify the fluid pressure by zone in the input file Syntax pressure value zone jxbegin jxend jybegin jyend jzbegin jzend fix or pressure value default Default None Explanation The pressure keyword uses an input format similar to other parameters specified by zone e g tortuosity permeability The keyword is follow
50. arl I Steefel 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 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 is assumed to be just uC the same as the advective flux at the ext 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 concentration boundary condition At the 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 zer
51. ase of an injection well Infiltration Keyword followed by a value s giving the infiltration or recharge rate the name of the geochemical condition to be used and the coordinate direction X Y or Z Syntax infiltration rate label coordinate Default None Explanation This keyword specifies an infiltration rate which is distributed over the boundary of the system at JX JY or JZ 0 The geochemical condition to be used is given by its label 78 CrunchFlow Manual Carl I Steefel Read_VelocityFile Keyword followed by a name for a file containing velocity vectors over the entire spatial domain Syntax read_VelocityFile filename format Backwards compatible read_velocity Default None Explanation This keyword provides a filename prefix for ASCII files containing velocity vector values over the entire spatial domain The extensions vx vy and vz are assumed to indicate the files containing the X Y and Z velocities respectively For example the syntax read VelocityFile copperleach SingleColumn would open and attempt to read the files copperleach vx copperleach vy copperleach vz using the SingleColumn format see Section 7 1 2 Units are set with the space units and time units keywords otherwise default units are m m yr The format of the file depends on the format chosen In the absence of a specification of the file format a single column SingleColumn format is assumed Velocities a
52. at Acad Sci 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 Kirkner D J and Reeves H 1988 Multicomponent mass transport with homogeneous and heterogeneous chemical reactions Effect of chemistry on the choice of numerical algorithm I Theory Water Resources Res 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 A C 1981 Rate laws in chemical reactions In Kinetics of Geochemical Processes ed A C Lasaga and R J Kirkpatrick Rev Mineral 135 169 Lasaga A C 1984 Chemical kinetics of water rock interactions J Geophys Res 4009 4025 Lichtner P C 1985 Continuum model for simultaneous chemical reactions and mass transport in hydrothermal systems Geochim Cosmochim Acta 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 143 165 Lichtner P C 1992 Time space continuum description of fluid rock interaction in permeable media J Geophys Res 3135 3155 Marsily G de 1986 Quantitative Hydrogeology Acad Press New York 440 p Patankar S V
53. ates is followed by a value for the logarithm of the mineral reaction rate that will be used for all of the minerals given in the MINERAL keyword block In this mode the code will not require that a kinetic rate law be explicitly given in the database for each of the minerals considered This option may be useful where the interest is in simply tracking a reaction path without detailed kinetic information 36 CrunchFlow Manual Carl I Steefel Gimrt Keyword followed by true or false which selects coupling method for reaction and transport Syntax gimrt logical logical is a standard FORTRAN logical true or false Default true Explanation The keyword parameter gimrt is followed by true or false or yes or no and is used to select a global implicit coupling of reaction and transport GIMRT if true or time splitting of reaction and transport OS3D if false Other restrictions on the coupling method may apply for example a global implicit method of coupling reaction and transport is currently allowed only up to 2 dimensions The selection of a global implicit method gimrt true eliminates any strict requirement that a CFL criteria be observed although other practical limitations on the time step either in terms of numerical stability or of time truncation error may apply See Steefel and MacQuarrie 1996 for a fuller discussion of the pros and cons of the two coupling methods Gimrt_pc Keyword followed by an option either bjaco
54. ay 9D 2 where m is the cementation exponent Dw is the diffusion coefficient in water and Dey is the effective diffusion coefficient in porous medium As defined here m 1 would be the formulation in porous media in the case where there is no tortuosity effect and the diffusion coefficient is corrected only for the cross sectional area actually made up of pores According to this formulation the formation factor would be defined as F Formation_factor Keyword followed by a value for the formation factor Syntax formation_factor value Default 0 Explanation This parameter provides an a formation factor which is defined as where F is the formation factor D is the diffusion coefficient in water and Dey is the effective diffusion coefficient in porous medium According to the definition of the formation factor in use in the marine geochemistry literature for example the cementation exponent should be set to 0 i e the porosity dependence of the effective diffusion coefficient is incorporated into the formation factor 73 CrunchFlow Manual Carl I Steefel Diffusion_activation Keyword followed by a value for the molecular diffusion activation energy in kcal mole Syntax diffusion_threshold value Default 5 0 default units of kcal mole Explanation This parameter sets the activation energy to be used in calculating diffusion coefficients at temperatures other than 25 C D_25 Keyword specifying a diffusion coe
55. bi or ilu Syntax gimrt_pe option Default bjacobi Explanation This sets the preconditioner method for the GIMRT option with the choices being a block Jacobi method bjacobi or ILU See the PETSc web pages at www anl gov petsc for further explanation Gimrt_pclevel Keyword followed by an integer value giving the level of fill to be used in preconditioning the global implicit reactive transport matrix Syntax gimrt_pclevel integer value Default 2 Explanation This sets the preconditioner level of fill for the GIMRT option as an integer value Using higher levels of fill than 2 may slow the speed of execution See the PETSc web pages at www anl gov petsc for further explanation Gimrt_solver Keyword followed by an option either gmres or bcgs Syntax gimrt_solver option Default gmres Explanation This sets the solver method for the GIMRT option with the choices being the GMRES method gmres or a stablized biconjugate gradient bcgs See the PETSc web pages at www anl gov petsc for further explanation 37 CrunchFlow Manual Carl I Steefel Graphics Keyword followed by an option either kaleidagraph tecplot or xmgr Syntax graphics option Default kaleidagraph for 1D tecplot for 2 3D Explanation This sets the style of output for the various graphics packages with Kaleidagraph being the default for 1D Tecplot the default for 2D The xmgr option can be used for either the X windows based xmgr or for gnuplot
56. c pressure gradient in the Y direction with depth increasing with increasing JY one could use pressure 1000 0 default pressure 0 0 zon 1 20 0 0 1 1 fix gravity 90 0 0 0 90 0 up With this input the code would time step until the hydrostatic pressure gradient was achieved in the absence of other forcings To initialize the system to hydrostatic pressure before time stepping begins use pressure 1000 0 default pressure 0 0 zon 1 20 0 0 1 1 fix gravity 90 0 0 0 90 0 up initialize hydrostatic true This will create initially hydrostatic conditions even if pressure boundary conditions or source sink terms are such that hydrostatic conditions will not prevail once time stepping begins 84 CrunchFlow Manual Carl I Steefel Initialize_hydrostatic Keyword followed by true or false which selects whether or not to initialize the domain to hydrostatic conditions Syntax initialize hydrostatic logical logical is a standard Fortran logical true or false Default False Explanation The keyword parameter initialize_hydrostatic is followed by true or false or yes or no and is used to select when true the option to initialize the pressure field to hydrostatic conditions before time stepping begins When set to true hydrostatic conditions will be initialized no matter what other pressure boundary conditions or source sink terms take effect once time stepping begins However a coordinate direction for the gravity vector must be
57. cay chain Pu gt Am gt Np can be represented by the following entries in the database Pu 241 decay 1 irreversible 2 1 0 Put 241 1 0 Am 4 83E 2 1 tot_Pu 241 1 0 Am decay 1 irreversible 2 1 0 Am 1 0 Np 1 60E 3 1 tot_Am 1 0 Since the reactions are designated as irreversible no equilibrium constants are needed Entries in Input File Aqueous kinetic reactions are entered in the input file using the reaction label that in each case begins a database entry So for the example above the AQUEOUS_KINETICS keyword block would look like the following AQUEOUS KINETICS Am _decay Pu 241 decay END No other options are currently allowed in the input file ION_EXCHANGE An ion exchange reaction can be described via a mass action expression with an associated equilibrium constant Vanselow 1932 Sposito 1981 Appelo and Postma 1993 The exchange reaction can be written in generic form as vACI aq uBX s lt gt uBCl aq vAX s where X refers to the exchange site occupied by the cations A and B The equilibrium constant Keq for this reaction can be written as Vanselow 1932 _ BCL AX ACL BX where the parentheses refer to the thermodynamic activities Several activity conventions are in wide use One possibility is the Gaines Thomas activity convention which assumes a reaction stoichiometry of
58. ceed by substituting Equation 9 into Equation 11 to obtain an expression for the total soluble concentration in terms of the primary species alone Substituting this result into Equation 33 gives essentially a kinetic version of the formulation A of Kirkner and Reeves 1988 N LIC ar KP 0C jal n X deut log Q K A k Tec K 1 0 m 1 j l No 04 C 14 CrunchFlow Manual Carl I Steefel where the y s 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 L into an algebraic expression and is discussed in the next section Operator Splitting 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 Cro n C At 2 followed by a reaction step Oee cj pee E At L C k 1 N N 33 R J 1 N G4 which is in turn followed by another 1 2 transport step Cr u m At 2 L C k 1 N N 35 In this scheme C are the concentrati
59. chFlow currently assumes in all cases that the reaction involves the dissolution of one mole of the mineral or solid phase in question that is it s stoichiometric coefficient is 1 Solid phase reactions specified in an input file involve two separate entries in the database file 1 a thermodynamic entry which gives the stoichiometry of the reaction the equilibrium constants as a function of temperature the molar volume of the solid phase and its molecular weight and 2 a kinetic entry which gives the rate law and rate coefficients for the reaction All mineral or solid aqueous phase reactions are considered to be kinetic in CrunchFlow except in the initialization procedure where equilibrium with a solid phase may be selected In practice equilibrium with respect to solid phases is attained by choosing rate coefficients that are large relative to transport rates Database Formats Thermodynamic Database Entries The thermodynamic database for solid liquid phase reactions begins with a separate line for each entry with the mineral or solid name being enclosed in single quotes This is followed by the molar volume of the mineral in units of cm3 mole 51 CrunchFlow Manual Carl I Steefel lt MineralName gt lt Molar Volume gt lt Number of species in reaction gt lt Stoichiometric coefficient lt SpeciesName gt lt Stoichiometric coefficient lt SpeciesName gt lt Log K array gt lt Molecular Weight gt
60. cks which may occur more than once and are used to set up boundary and initial conditions and source terms In each case they provide information on the actual concentrations of primary species and minerals needed to carry out a distribution of species calculation Only after a particular geochemical condition is read and processed is it allocated as a particular boundary or initial condition or as a source term Geochemical conditions are input via a CONDITION keyword block The CONDITION string is followed in each case by a label which identifies the condition elsewhere in the input file Most of the information in a condition block pertains to aqueous and solid phase concentrations however certain parameters may be input which apply to the condition as a whole These are discussed in the following section Units Concentration units may be set using the units keyword within a particular condition The change in concentration units only applies to the condition within which it is set The format for input of different units is units concentration units The default concentration units for aqueous species is molality moles solute per kg water The additional options include Concentration Units Keyword Designator moles kg water mol kg millimoles kg water mmol kg micromoles kg water umol kg micromol kg parts per million ppm Table 3 Options for specifying concentrations in geochemical conditions
61. d use a timestep smaller than a Courant number of approximately 1 2 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 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 heterogeneous 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 differential operator L U 2o Ve u pv 32 if we drop the density and mass fraction of H O terms for simplicity and write the governing differential equations in terms of the total component concentrations as L U R J 1 N 33 min where the reaction term R 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 pro
62. d only in GIMRT mode and flux In the case of a flux specification the code uses the value and direction of the advective flux that is provided if the flux is into the domain then the upstream ghost cell will be used if out of the domain then there is no effect For a pure diffusion problem 1 e with no advective flux via flow across the boundary of the domain a flux specification will result in a no flux boundary In the case ofa Dirichlet boundary in GIMRT mode this will cause both advective and diffusive flux to be calculated To use a Dirichlet condition in OS3D mode requires that an initial condition be fixed just inside the boundary of the problem TRANSPORT The TRANSPORT keyword block reads in the inputs controlling transport of species Flow is treated in the FLOW keyword block The most important inputs concern molecular diffusion and hydrodynamic dispersion along with transport properties of the medium e g the cementation factor or tortuosity that affect these processes Fix diffusion Keyword followed by a value for the molecular diffusion coefficient in the space and time units set by the space units and time units keywords within the TRANSPORT block Syntax fix diffusion value Default 0 0 default units of m yr Explanation If this parameter is set then additional instructions to calculate diffusion calculate diffusion based on a 25 C value and an activation energy are ignored If no species specific dif
63. directly into the MINERALS keyword block However unlike similar operations for secondary species and gases the minerals specified in the input file also require kinetic information which resides in a separate part of the database from the thermodynamic entries for minerals or solid phases Since the kinetic database is not as complete as the thermodynamic database letters written to funding programs within DOE lamenting this situation will be appreciated some minerals found during the sweep may not have kinetic database entries yet This will be remedied in future releases by 1 augmenting the kinetic database for CrunchFlow and 2 by sweeping the kinetic database in addition to the thermodynamic database for solid phase reactions which can be loaded directly as kinetic reactions The database sweep will tend to find a relatively large number of phases which may or may not actually form at the temperature of interest It is common for example to find various forms of garnet supersaturated at 25C geochemists and petrologists know that these 56 CrunchFlow Manual Carl I Steefel phases precipitate only at high temperature Thus editing of the mineral list generated by a database sweep is essential Without any additional information the code will assume that the default reaction is to be loaded if only the name of the solid phase is given that is it will search for a kinetic entry of that name labeled as default Following
64. domain Syntax read TortuosityFile filename format Backwards compatibility read tortuosity Default None Explanation This keyword provides a filename for a file containing tortuosity values over the entire spatial domain The tortuosity is considered here to be a number less than 1 some formulations divide by the tortuosity so that it is always greater than 1 The tortuosity 1 as used in CrunchFlow therefore can be viewed as a parameter which allows an effective diffusion coefficient in porous media Dpy to be calculated from the diffusion coefficient in pure water D and the porosity according to Dy 7D The format of the file depends on the format specified after the file name but the values of the tortuosity are assumed to be listed with NX varying first then NY and then NZ For those familiar with FORTRAN and using the SingleColumn format the default the actual source code is Do jz 1 Do jy tortuosity jx jy jz Tortuosity Keyword used to specify tortuosity by zone in the input file 75 CrunchFlow Manual Carl I Steefel Syntax tortuosity value zone jxbegin jxend jybegin jyend jzbegin jzend or tortuosity value default Default 7 0 Explanation The tortuosity keyword uses an input format similar to other parameters specified by zone e g pressure permeability The keyword is followed by the value of the tortuosity which is then followed by either the label zone or default
65. e not neglected Future versions of the code will have the capability of reading all of the potentially relevant species from the database The only 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 A1 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 O among the list of gases if O is to be present in the system since much of the EQ3 database is written in terms of O 2 aq Tf 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 H O presents a special case since in many cases it is present in such abundance that its activity can be taken as unity If one wants the activity of H O to be 1 then it should be left out of the list of primary and secondary species The code will still balance the r
66. e ee eee eee ee eese sesso sesso seesssseesesesesesess 70 Exchanger Concentrations scsec vesessssecssesecsovdvesvsevnssvncnsesvonsssnesvssnstvsvntesveenevesansvevesteenesosveecs 70 INITIAL CONDITIONS i ccantpiccsiveitst on einai ead amieuenneenmwen TL BOUNDARY CONDITIONS 4i cssscsceszsrecasthcieebioteccksneerecnti austaelainane tte cee 12 Fix diffusion nine D Calculate diffusion TEE S Cementation_exponent ssssssecssecesocesocesoocssocessocesocescosesocssocessocesocssoosssocsssecesocesoossssesssese TO Formation factor eue ea ORI Gta md tetuer RON UA Diffusion ACC VAG ON mL S 0000000000000000000000000000000000000000000000000 0000000000000000000000000000000000000000000000000 Gas LITER SITE eoe order aeo be dee n ee EOE AYER CUN Aer dSE EOD nete ti Rare diis eite LA Read Vlil indo Pep M 2D FOPODOSIE ASSE ATE E E FO Anisotropy PATO Y cdcdccidedasiscestedesiuteviedecddacesdetsseedesdsddesededeaisescvicdacseteasstdonddesdseuaiatiassecesice dO PNE eet tyr LI COP P cess veassosssudenveiiescdadsveasdesecudenseiescessnstvsessesetsssefseedsodsvens 10 Threshold porosity ssssscissssscsessssssssosassssssssssssssssssssesssossssisossssessovss sssssssssssssssisssossssossssss ss d J Tortuosity b6 0Wuinidenetkssace kai tke aede awchb suben d aiie soasa s oiae deo Ceca b sasons as issant kasseta
67. e linked to properties like the porosity which in turn may be affected by mineral dissolution and precipitation Interface permeabilities are calculated as harmonic means between adjacent cells For example between nodes i and i the harmonic mean of the permeabililty is given by 2 kk i i l FS k ki l Since permeabilities need to be defined at the boundaries to determine the amount of flow across them ghost cell X permeabilities are required at the 0 and nx 1 nodes and so on For 1D problems permy and permz files are not needed For those familiar with FORTRAN the actual source code in the case where a SingleColumn format is specified is given by compare to the ContinuousRead format in the case of the read_VelocityFile keyword DO jz 1 nz DO jy DO jx DO jz 1 permy jx jy jz DO jz 0 Tx n 1 ny 54 permz jx jy JjZ 81 CrunchFlow Manual END DO END DO END DO Permeability_x Carl I Steefel Keyword used to specify the permeability in the X coordinate direction by zone in the input file Syntax permeability_x value zone jxbegin jxend jybegin jyend jzbegin jzend or permeability_x value default Default None Explanation The permeability_x keyword uses an input format similar to other parameters specified by zone e g pressure tortuosity The keyword is followed by the value of the X permeability which is then fol
68. e master species specified by the keyword parameter master This parameter can be used to either control the time truncation error or to limit the increase in time steps for the sake of numerical stability With a time tolerance set very high the code will normally increase the time step by a factor of 2 every time step until the maximum time step is reached 43 CrunchFlow Manual Carl I Steefel OUTPUT The keyword block OUTPUT controls when and what will be written to disk If the OUTPUT block is not given then no time stepping will occur see keyword parameter spatial_profile The keyword parameters included in the OUTPUT keyword block are time units spatial profile time series time series print time series interval MakeMovie The time units keyword parameter has the same definition and role as in other keyword blocks where it appears It specifies the time units for such keyword parameters as spatial profile and it determines the units used in the various output files Example OUTPUT time units hours spatial profile 250 time series 50 time series 100 time series print Cs pH time series interval END In this example a single spatial profile field variables as a function of space is written at 250 hours The code will time step until 250 hours is reached and then stop A total concentration of cesium will be tracked at JX 50 and JX 100 nodes 50 and 100 in the X direction Output to the time series
69. eactions properly as if H O were present If on the other hand one wants to carry out a mass balance on H O then it should be included in the list of primary species and it will be treated like any other 27 CrunchFlow Manual Carl I Steefel species except for the treatment of its activity coefficient GIMRT assumes that the activity of H O is given by its mole fraction in the solution Redox Reactions There are many possibilities for choices of primary and secondary species to represent redox reactions GIMRT writes 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 There are cases where electrons need to be included as in corrosion problems where an actual electron flux occurs When the solution is electrically neutral however no special provision for redox species and reactions is necessary when whole cell reactions are used in contrast to statements by Yeh and Tripathi 1989 although the inclusion of redox reactions implies that we have included an additional balance equation representing the conservation of exchangeable charge Since O included in the list of gases in the input file if redox reactions are to be considered A gas however can
70. ecies concentrations and the equilibrium constant for the reaction using a mass action expression Thus the species can be mathematically eliminated as an independent unknown from the conservation equations In contrast where the reaction is kinetically controlled a differential as opposed to algebraic equation is needed to relate the species to the other species in the problem Database Format As with the primary species described above a distinction must be made between secondary species in the database versus secondary species specified by the user in the input file The database secondary species block follows the primary species block and is initiated with the line marking the end of the database primary species block End of primary 0 0 0 0 0 0 and terminated by End of secondary 1 0 0 Ou Oi Ow Q0 0 0 0 0 0 0 0 Database secondary species are written in every case as the destruction one mole of the secondary species Using the standard convention therefore the stoichiometric coefficient for the secondary species is assumed to be 1 even though this is not explicitly given in the database The format for a database secondary species is SpeciesName gt Number of species in reaction integer gt Stoichiometric coefficient gt lt SpeciesName gt lt Stoichiometric coefficient lt SpeciesName gt lt Log K array Debye Huckel size parameter Charge Molecular we
71. ecified within the RUNTIME keyword block OvershootTolerance Pc Keyword followed by a real number giving the tolerance for overshooting a mineral volume fraction Syntax OvershootTolerance value value is a positive real number Default 1 0E 05 Explanation This keyword sets the tolerance for overshooting a zero volume fraction This is necessary because of the way in which CrunchFlow handles mineral abundances For any one time step mineral abundances and surface areas are assumed fixed qmineral volume fractions therefore are updated with an explicit scheme at the end of the time step This approach however makes it possible to overshoot the available mineral concentration which can result in overall mass balance errors The tolerance set here specifies the maximum overshoot in units of m mineral m porous medium If the tolerance is exceeded the time increment delt is cut and the time step is repeated Keyword followed by an option either ilu jacobi or direct Syntax pc option Default ilu Explanation This sets the preconditioner method for PETSc to solve the flow or diffusion equations It does not apply to the solution of the global implicit reactive transport equations in the GIMRT option The choices include i u with a partial fill set by the keyword pclevel jacobi and direct The direct option carries out a direct solution of 39 CrunchFlow Manual Carl I Steefel the matrix without iteration 1 e it is
72. ed by the value of the fluid pressure in units of Pascals which is then followed by either the label zone or default If zone is specified then the code expects the X Y and Z coordinates in the form of a beginning JX and ending JX a beginning JY and an ending JY and a beginning JZ and an ending JZ in each case separated by a hyphen If default is specified instead of zone the code will initialize the entire spatial domain 1 e all grid cells to the value More than one specification of pressure can be and normally is provided for example pressure 1000 0 default pressure 0 0 zone 1 20 21 21 1 1 fix Zones specified later in the sequence overwrite zones specified earlier in the example above the default specification sets the fluid pressure to 1000 Pa initially in all of the grid cells while the next keyword sets the pressure at JY 21 to zero In the case of a 20 by 20 X Y simulation the grid cell 21 corresponds to node NY 1 and is therefore a ghost cell This is the normal way in which fixed pressure boundary conditions are set Following the zone coordinates with fix instructs the code to lock the value in for the course of the simulation not including fix would mean that the pressure is only set as an initial condition after which it is free to evolve according to the physics of the problem The convention adopted here is that all three dimensions must be specified even in a one dimensional problem To create a hydrostati
73. ees OU Entries TIE Eileen dO Re entente a eee O2 SURFACE COMPLEXATION ieitiickteseo erben due Reeos reb u baba pte eto cosa qiu rer antecctidestecns OS GEOCHEMICAL CONDITIONS Equilibrate surfaCC s sssissssessssssssssisssressssses ssossossssseisorsssssassssesussstoozs sssseorsssvssozisse sekss ss OF 00000000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000 PKAN AATE PEETA Ma ER cH ERR NUR Udo Eel pe qa EIN DERE abo e TRE ADT U AREA AEREA OS Set DOPOSIUS aei er p eder T r tM E etn E sni cte ters Pda OD SEE Satu FIIO cierre seco debe bi i vEUs ede bed seul ka vba doses sdotdevidcsadecessecsedkstbavdeiattoussdosdssesisee OD Types of Constraints Aqueous Species eee ee eee creen eee eee eee ee eese ense ee ssseessssees D ssessssceeee 06 scl Species Activity Constraint i ee eetieeto tovto tatit Se proe pe PUER PE MEE Erie ori Mi ecieae tese eiicieiee OT Gas CODSEPAID I ck rores I ERIS eR Creo EPUM Fee EFE R FEES Ser E PV PER AR bore ot RERO E PRSE Po reses OF Mineral Constraint nr Charge C OMS UAE Sec OT Solid Phase Inputs soa casssccceesacdenceesxsaae BRAIN I I ASURRENEINNEIN NU cagnetuaneacsesndoeeaseanmaems prd erre OS Mineral Volume Fractions sisse 00 Mineral Surface Area Options METUO EU IN DUO MEE NOE ES S esieeeeeese OD Surface Hydroxyl Concentrations ecce ee
74. efel NZ For those familiar with FORTRAN and using the SingleColumn format the default the actual source code is Do jx 0 nx READ 52 FluidBuryX jx jy jz SolidBuryX jx jy j2z End Do 89 CrunchFlow Manual Carl I Steefel Acknowledgements Many people have delayed the completion of the manual When I think of anybody who has contributed to speeding its completion I will list them here Literature Cited 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 237 285 Aris R 1965 Prolegomena to the rational analysis of systems of chemical reactions Arch Rational Mech and Anal 81 99 Bear J 1972 Dynamics of Fluids in Porous Media American Elsevier 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 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 465 480 Daus A D and Frind E O 1985 An alternating direction Galerkin technique for simulation of contaminant transport in complex groundwater syste
75. eft like this Fe and Fe would be treated as if they were separate elements with no reaction relationship linking their concentrations This is the so called redox disequilibrium option found in a number of the reaction path codes In CrunchFlow however one can add a kinetic reaction linking these species using the AQUEOUS KINETICS keyword block 59 CrunchFlow Manual Carl I Steefel Database Formats Unfortunately the aqueous kinetics inputs have not yet been shifted over to the more readable format of the mineral kinetics section In addition the aqueous kinetic reactions do not currently make use of the equilibrium constants and reaction stoichiometries found in the thermodynamic portions of the database This is nominally because kinetic pathways may have a specific stoichiometry which differs from that found in the thermodynamic entries As in the case of mineral water reactions parallel reaction pathways may be operating each associated with its own different reaction stoichiometry In the initialization speciation procedure since equilibrium is assumed only a single reaction stoichiometry is needed Steefel 2000 Any reaction stoichiometry in terms of the primary species can be used and will give identical results mathematically The aqueous kinetics section follows the MINERALS section of the database and is begun with the keyword phrase Begin aqueous kinetics and is terminated with the keyword phrase End of aque
76. ely for the bulk and specific surface area options If the surface area option is left off the code assumes by default that the trailing numerical value refers to the bulk surface area Specification of a bulk surface area will lead to a calculation of specific surface area according to Apuk V m A specifie 2 MW where V is the molar volume of the solid phase is the initial volume fraction and MW is the molecular weight of the phase The specific surface area is then used together with other 68 CrunchFlow Manual Carl I Steefel surface complexation parameters to calculate concentrations of surface hydroxyl sites per volume porous medium The two surface area options lead to different methods for computing mineral surface area as a function of time Bulk surface area In the case where the bulk surface area is specified the reactive surface area of a solid phase as a function of time is calculated according to 2 3 2 3 Apuk ATE e orna VT dissolution 2 3 A 157729 precipitation where refers to the porosity and 4 refers to the individual mineral volume fraction The inclusion of a 2 3 dependence on porosity is chiefly to ensure that as the porosity goes to 0 so too does the mineral surface area available for reaction This formulation is used primarily for primary minerals that is minerals with initial volume fractions gt 0 For secondary minerals which precipitate the value of the initi
77. endence Rate laws for tst type reactions take the following form rate mol kgw yr gt Number of species dependences gt lt SpeciesName gt lt Exponent gt lt SpeciesName gt lt Exponent gt 60 CrunchFlow Manual Carl I Steefel As an example consider the kinetically controlled oxidation of Fe by molecular oxygen which uses two parallel reactions one pH dependent and one pH independent following the rate law proposed in Singer and Stumm 1970 FelI oxidation 2 tst 5 1 0 Fe 0 500 H20 0 25 O2 aq 1 0 Fet 1 0 H 8 4887 T 5b3e 06 3 Vtot Feet 1 0 02 aq 11 0 AHE 2 0 41 4848 2 tot Fett 1 0 O2 aq 1 0 As in the case of solid liquid phase reactions the reaction rate can be made to depend on the total concentration of a species Fe in this case by prepending the string tot to the species name Monod Monod reactions take a slightly different form than do TST type reactions The Monod option has a field for inputting various Monod terms that is the dependence of the reaction rate on electron acceptors and or electron donors R E K nax o C K half where C is the concentration of the electron acceptor or donor Krais the half saturation constant in moles kg water and kmax is the maximum rate of the reaction when the concentration of the electron donor or acceptor gt gt the half saturation constant The specification of Monod
78. eobacter and 2 the abiotic reduction by hydrogen sulfide FeOOH s 1 925H 0 033NH 0 208CH COO gt Fe 0 033C H 0 N 0 25HCO 1 6H O designated in the input file and database as Iron DIRB while the reaction 2FeOOH HS 5H 2Fe S 4H O could be designated as Iron H2S The code thinks these are two separate phases because of the different entries and stoichiometries unless the associate keyword is used If one mineral is associated with another then only the volume percentage and reactive surface of the associated mineral is updated So in the example above the following usage 58 CrunchFlow Manual Carl I Steefel MINERALS Iron DIRB Iron H2S associate Iron DIRB END would result in the dissolution according to the Iron H28 kinetic pathway affecting the volume fraction and surface area of the Iron _DIRB phase To see that this is working users can check to see which of these phases is being updated as a function of time in the volume out files AQUEOUS KINETICS The keyword block and database entries for aqueous phase homogeneous kinetics currently have a simple form which will be changed in later releases to conform more closely to the format used for mineral water or solid liquid reactions Aqueous phase reactions must be made up entirely of primary species the concept of a secondary species in CrunchFlow is restricted to reactions where an equilibrium relationship b
79. es a FERMER eda a Rose EIH Fo Das eu Daans CERTI VR Cra e E eL a tad ku basati oaei 2 L wrredounMinitsiloi em P G Input File Formats biis d SingleColumn Example ESEE P ContinuousRead Example 4 cxstevccsnciecencicagudacoesucoraceuvedecsesanseserassoendessengsvonsedeaneaterdovenaeses OL Fullbormat Example Soon EGO RN OUR IU Ota Ed ee Unform tted Dori ru rm NR OO 00000000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000 000000000000000000000000000000000000000000000000 0000000000000000000000000000000000000000000000000 aa Courant number e e eee ee esee esee eese esee eese ee sese ee sese se sese se sese sese se sess esses sess sese seseses sess cesse OF 0000000000000000000000000000000000000000000 Database 5 e atre Se er na EEE gena Pee eR e Prae tosses caseesceesessesecsecocdccsessesceseescoecesteaseeceason OD Database SWEEP onde nseteridiscieseua e evi dean cd otl oe nra da eid esaiok a dk Red V ska TEE FER R dE Pe eb AREE Pd a eed OO Debye Huh eT Density mod le 4 45 nudsidscibn deer kite oe ona Eesob avkh ess n reor ibid ede Do ania ernie OD Dissol tion MaX qood ect bit rn CE RRBRROEDEDEN US Xn Cod UE EE UTOR MARER C
80. es present in the secondary species The U s therefore 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 10 CrunchFlow Manual Carl I Steefel As noted above the choice of H O H 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 E 4 ce ERANS jt 1 0 OH Al OH so that the total concentrations now become Uso Cao 4 gus t Cu 24 ALT OH 23 U ys Cy 4C 7C 25 Al OH U Com C C6 Al OH B Al OH Writing the chemical reactions differently which mathematically is referred to as a change of basis therefore changes the total concentrations of H O and H even though the total concentration of the Al bearing species is the same Reaction Rate Laws We use a kinetic rate law based on the assumption that attachment and detachment of ions from mineral surfaces is the rate limiting step 1 e a surface reaction controlled rate law It does not mean however that one cannot obtain overall transport control on the mineral dissolution or precipitation rate since this depends on the magnitude of the reaction rate relative
81. etween the species in the reaction applies In an aqueous phase kinetic reaction species which might have been secondary species if equilibrium applied need to be moved to the list of primary species or the reaction will not be loaded Consider the case of the Fe ID Fe IIT redox couple In the case where the oxidation of Fe II by molecular oxygen occurs rapidly and therefore can be assumed to be at equilibrium one could use a list of primary species given by the following to describe the system PRIMARY SPECIES H 7 Fet O2 aq END followed by specification of Fe in the list of secondary species This would then load the reaction neglecting the trailing equilibrium constants Fet 4 0 50 H20 0 25 O2 aq 1 00 Fe 1 00 Ht along with its associated equilibrium constants Fe would then become a secondary species contributing to the total concentrations of all of the primary species in the reaction not only Fe but also O2 aq and H A fuller discussion of this can be found in Steefel and MacQuarrie 1996 or in Bethke 1996 However if one wanted to treat this reaction as a kinetic one then Fe needs to be listed as a primary species along with the others PRIMARY SPECIES H i F 02 aq F END If this is done then the reaction involving Fe given above will not be loaded from the thermodynamic portion of the database L
82. fficient at 25 C for a specific aqueous species Syntax D_25 species_name value Default value set in fix_diffusion keyword Explanation This keyword provides a diffusion coefficient in the units based on the space_units and time_units keywords set in the TRANSPORT block When even one species specific diffusion coefficient is given the code GIMRT only solves the Nernst Planck equation which includes an explicit calculation of electrochemical migration In this case any species not listed with this keyword are given the diffusion coefficient specified in the fix_diffusion keyword Gas_diffusion Keyword followed by a value for the gas diffusion coefficient Syntax gas_diffusion value Default 0 0 default units of m yr Explanation This keyword sets a value for the gas diffusion coefficient in the units given in the time_units and space_units keywords set within the TRANSPORT block Dispersivity Keyword followed by a value for the gas diffusion coefficient Syntax dispersivity longitudinal value lt transverse value gt Default 0 0 default units of m Explanation This keyword sets a value for the longitudinal first entry and transverse second value dispersivities Space units may be changed with the space units keyword set within the TRANSPORT block If the problem is one dimensional the transverse dispersivity may be omitted Constant tortuosity Keyword followed by a value for the tortuosity or by a string giv
83. 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 18 CrunchFlow Manual Carl I Steefel direction In the case of flow in the positive x direction 1 e u gt 0 the estimated concentration would be Cric C C B r Er 56 where 24r f r Max n Min 22r zu 57 and C C C C r i i l i l i 5 8 X t x x t x Similar calculations are performed for negative velocities and the method is generalized for multiple dimensions 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 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 numerical scheme second order accurate in space and time is applied to the diffusion equation The diffusion 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 fo
84. fusion coefficients are specified see D 25 keyword With setting of this parameter all aqueous species are given the same diffusion coefficient and only Fick s law is implemented in the calculation the electrochemical migration term is zero when all species have the same diffusion coefficient in any case If any species specific 12 CrunchFlow Manual Carl I Steefel diffusion coefficients are set then the full form of the Nernst Planck equation is solved see section XXX and the value given by fix diffusion is used for any species not given a specific diffusion coefficient value with the D_25 keyword Calculate_diffusion Keyword followed by a value for the molecular diffusion coefficient in the space and time units set by the space_units and time_units keywords within the TRANSPORT block Syntax calculate_diffusion value Default 0 0 default units of m yr Explanation This parameter is used to calculate a diffusion coefficient at 25 C Diffusion coefficients at other temperatures are calculated with the value set in the diffusion_activation keyword The fix_diffusion keyword over rides this keyword for clarity both should not normally be set Cementation_exponent Keyword followed by a value for the cementation exponent in Archie s Law Syntax cementation_exponent value Default 1 0 Explanation This parameter provides an exponent to be used in Archie s Law to specify a porosity dependent tortuosity according to D
85. gz jx jy 22 x 1 nx jy 1 ny jz 0 nz uj ui ug FullFormat Example DO jz 1 nz DO jy 1 ny DO jx 1 nx READ 52 xdum ydum zdum tortuosity jx jy jz END DO END DO END DO 32 CrunchFlow Manual Unformatted Example READ 52 qx TITLE Carl I Steefel The TITLE keyword block is used to provide a title for a particular simulation It is echoed to the standard output during a run and written to the output file crunch out The title will be read as a character string up to 264 characters long Example TITLE Lichtner s 2D copper leaching problem END RUNTIME The RUNTIME keyword block contains a number of parameters for running a particular simulation All of them are optional as is the RUNTIME keyword block itself since defaults values will be used if none are provided by the user All of the parameters with the exception of the species name given in the keyword parameter master are case insensitive Within the RUNTIME keyword block the following keyword parameters are possible coordinate correction max courant number database database sweep debye huckel density module dissolution max fix saturation generic rates gimrt gimrt pc gimrt pclevel gimrt_solver graphics hindmarsh lag_activity later inputfiles master pc pc level precipitation max reaction path restart save restart screen output solver speciate only steady state timestep max ti
86. h these names for a particular reaction to be loaded This name must also match an entry in the thermodynamic part of the database although this requirement may be relaxed in the future for the special cases of irreversible reactions where thermodynamics may be irrelevant The code also allows for multiple parallel reactions multiple reactions running in parallel which affect the same solid phase through the database entry on the line following the solid phase name This line begins with the string label and is to be followed with the label of the rate law The following third line begins with the string type which may include currently up to five different types of rate law see below The fourth line begins with rate 25C The fifth line gives the activation energy for the reaction While there are different options after the fifth line depending on the type of reaction involved see below all of the solid liquid phase kinetic entries have a common input format for the first five lines following the separator Calcite label lt name gt type lt type of reaction gt rate 25C lt rate at 25C in log base 10 units gt activation lt activation energy kcal mole gt The label option can be used to construct an overall rate of reaction that is the sum of a number of parallel reaction pathways An example would to use several database entries for Calcite 52 CrunchFlow Manual Carl I
87. he GASES keyword block specifies the user s choice of gases for a particular problem In the case of fully saturated flow gases are only used in the initialization procedure at the present time since their mass is not explicitly accounted for Future releases of the code will have an option to specify a kinetic reaction with a gas phase which could be used for example to continually equilibrate a solution with the atmosphere Clever users will recognize that the same effect can be obtained presently by specifying a fictional mineral in the MINERAL keyword block whose equilibrium constant corresponds to the aqueous concentration at which a species would be in equilibrium with a particular partial pressure of the gas phase an example is given below During the initialization of the various geochemical conditions the option is available for equilibrating the solution with one or more gases at specified partial pressures In previous 50 CrunchFlow Manual Carl I Steefel versions of the code it was necessary to input the gas O2 g in order to get most of the redox reactions to load properly In the current version however the specification of O2 g is optional since the code will automatically load it if the aqueous species O2 aq is found in the primary or secondary species list For fully saturated problems therefore gases need to be input only if they are used in constraining a geochemical condition For unsaturated transport and gas
88. he name of the exchanger must correspond to one of those identified in the ION EXCHANGE keyword block Solid density may be specified one of three ways 1 Direct specification of the bulk solid density kg m SolidDensity Value 70 CrunchFlow Manual Carl I Steefel 2 Specification of the solid solution ratio g kg water g L SolidDensity ss_ ratio Value 3 Calculation based on the sum of the mineral volume fractions SolidDensity CalculateFromMinerals Exchange on a specfic mineral Here the total number of equivalents of exchange sites is calculated from the combination of the mineral volume fraction and the cation exchange capacity per gram of that mineral The solid density therefore is not relevant the keyword is therefore ignored although the values of the porosity and liquid saturation affect the result see below Mineral densities are calculated from the combination of molar volumes and molecular weights in the database As in the case of exchange on the bulk material the CEC is specified with a similar format but here the units are equivalents per gram mineral Exchanger cec Value INITIAL CONDITIONS In order to run a reactive transport or a reaction path calculations initial conditions must be set Even in the case where the desired result is a steady state simulation in which initial conditions are irrelevant the code must still achieve this steady state by stepping through time Only with the
89. ies witha O2 aq O2 g 0 20 gas Mineral Equilibrate primary species witha O2 aq Pyrite mineral Charge Charge balance the solution using Na charge the primary species Table 4 Constraints available for specifying aqueous species concentrations The most commonly used constraint is the total concentration constraint Analyses of natural waters generally report total concentrations rather than individual species concentrations Or someone may know how much total mass of a chemical element they added to an experimental system While specification of individual species concentrations is relatively uncommon specification of species activity is common especially for hydrogen ion pH and for molecular oxygen O2 aq Total Concentration Constraint The total concentration constraint signals to CrunchFlow to carry out a total mole balance calculation If the balance is over the aqueous phase only the default the equation takes the form Ns T C 4C where Tj is the total concentration of the primary species Cj is the individual primary species concentration v j is the stoichiometric coefficient and Ci is the concentration of the Ns secondary species If the keyword equilibrate_surface is specified either as a global parameter within a condition or for an individual primary species then the mole balance equation also includes the mass of the primary species incorporated in ion exchange sites and as surface c
90. ight gt In this format the number of pairs of stoichiometric coefficients and species names is determined by the preceding value in Number of species in reaction if there is a mismatch a read error will result The length of the Log K array is given by the number following temperature points in the first line of the database file Normally eight temperature points are used following the EQ3 database format for example temperature points 8 0 25 60 100 150 200 250 300 In this example log K values for the reaction will be provided at eight temperatures between 0C and 300C 49 CrunchFlow Manual Carl I Steefel An example of an entry for the species CO2 aq is given by CO2 aq 3 1 00 H20 1 00 H 1 00 HCO3 6 5804 6 3447 652684 56 3882 6 7235 7 1969 7 7868 8 5280 3 0 0 20 44 0098 Input File Entry of Secondary Species The secondary species list provides the additional aqueous species used in a reactive transport problem If the user selects the database sweep option in the RUNTIME keyword block however the code will automatically load all possible secondary species ignoring the user s list However when this option is selected the code will stop after speciating the geochemical conditions To use the set of secondary species generated by the database sweep option the user needs to formally copy and paste the secondary species list into the SECONDARY SPECIES keyword block and tu
91. ile are given in the Output Variables section below Time_series Option to specify nodal location of a time series Syntax time series filename JX JY JZ for a 1D problem time series filenamel JX JY JZ for a 2D problem time series filename3 JX JY JZ for a 3D problem where filename is the file name to be used for the time series and JX JY and JZ are integers giving the node number for which the time series will be tracked Default There are no default values Explanation This keyword parameter specifies that a time series should be recorded at the node specified Multiple specifications of time series may occur each one giving a different filename and node number For one dimensional problems the Y and Z coordinates are optional If in any case a node location is given that is greater than the coordinates specified in the DISCRETIZATION block execution will terminate and an error statement will be written e g JX gt NX the number of nodes in the X direction Time series print Option to write out a time series for one or more solutes Syntax time series print speciesI species2 species3 amp speciesn where speciesn is a character string corresponding to an aqueous species name Default If fime series print is not specified no tracking of solutes will be done If time series print is not followed by a species name a default of all will be assumed Explanation This keyword parameter gives the species to be t
92. ing the name of the tortuosity option Syntax constant tortuosity value or constant tortuosity option Default 0 74 CrunchFlow Manual Carl I Steefel Explanation This keyword sets a constant value for the tortuosity which is then applied throughout the spatial domain Definitions of the tortuosity vary in most cases it does not include the porosity that multiplies the diffusion coefficient and this is the convention adopted here The tortuosity is also considered here to be a number less than 1 some formulations divide by the tortuosity so that it is always greater than 1 The tortuosity 1 as used in CrunchFlow therefore can be viewed as a parameter which allows an effective diffusion coefficient in porous media Dpy to be calculated from the diffusion coefficient in pure water D and the porosity according to Dy 7D If an option rather than a numerical value is given the code checks the option specified against the list currently implemented At the time of this manual the only option that is accepted is CostaRica which is a porosity tortuosity option to published soon With an option set the tortuosity is calculated dynamically at run time so the keyword constant tortuosity is a bit of a misnomer here but it avoids a read of tortuosity from an input file read tortuosity or by field tortuosity Read TortuosityFile Keyword followed by a file name containing tortuosity values over the entire spatial
93. ion of a mineral 0 its surface area also 0 In addition for both dissolving and precipitating minerals the term g A requires that the surface area of a mineral in contact with fluid gt 0 when the porosity of the medium gt 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 B dissolution 2 3 H precipitation h which implies that if a secondary were to begin to dissolve its surface area would equal its initial surface area A When V 1 if the porosity were constant An Ano 62 Program Structure 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 series 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 N xM algebraic equations in order to obtain the corrections to the component concentrations 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
94. ions of some of the species will not change from one representation to the other for example total HCO will equal total CO aq However other total concentrations will change including that of H Use of the Fe Fe redox couple is an alternative way of handling redox than is given by the use of Fe and O aq Note that the definition of total elemental Fe will change depending upon which formulation is used Example 1 oe PRIMARY SPECIES d Nad K Example 2 E PRIMARY SPECIES Nat Al t K F Ca SiO2 aq H HCO3 A S04 p 48 CrunchFlow Manual Carl I Steefel Fe Cis Si02 aq Fe CO2 aq END SO4 SECONDARY SPECIES The SECONDARY SPECIES keyword block gives the user s choice of secondary species for a particular problem Secondary or non component species are those for which an equilibrium reaction relationship is assumed with the primary species in the problem In other words the secondary species can be written in terms of the set of primary species While the secondary species could be viewed as any species for which there is a reaction relationship with the primary species whether the reaction is at equilibrium or kinetically controlled in this code the secondary species are used to imply on species where the reaction is assumed to be at equilibrium When a particular reaction is at equilibrium a secondary species can be written algebraically in terms of the primary sp
95. irection relative to the value in the X direction Syntax Anisotropy_ratioY value Default 0 Explanation This keyword defines the ratio of the diffusion coefficient in the Y direction relative to the diffusion coefficient in the X direction The anisotropy term therefore is dimensionless If the anisotropy factor in the Y direction is defined as y then the effective diffusion coefficient in the Y direction JY is given by Dy 91D Anisotropy ratioZ Keyword followed by a value defining the anisotropy ratio in the Z direction relative to the value in the X direction Syntax anisotropy ratioZ value 76 CrunchFlow Manual Carl I Steefel Default 7 0 Explanation This keyword defines the ratio of the diffusion coefficient in the Z direction relative to the diffusion coefficient in the X direction The anisotropy term therefore is dimensionless If the anisotropy factor in the Z direction is defined as z then the effective diffusion coefficient in the Z direction JZ is given by Dox 9tD 7 Threshold_porosity Keyword followed by a value defining a porosity threshold separating two different tortuosity values Syntax threshold_porosity value Default 0 0 Explanation This keyword sets a threshold porosity either side of which the tortuosity takes on differing values given by the keywords tortuosity_below and tortuosity_above Tortuosity_below Keyword followed by a value giving the tortuosit
96. le co2stream ContinuousRead would open and attempt to read the files co2stream vx co2stream vy co2stream vz using the ContinuousRead format see Section 7 1 2 Units are set with the space units and time units keywords otherwise default units are m m yr The format of the file depends on the format chosen In the absence of a specification of the file format a single column SingleColumn format is assumed Velocities are defined at grid interfaces not at the center of the nodes so NX 1 X velocities are needed in the X coordinate direction NY 1 Y velocities in the Y coordinate direction and NZ 1 Z velocities in the Z coordinate direction CrunchFlow uses the convention that velocities for a particular grid cell are defined at the right hand interface the so called right hand rule So for example the gas velocity qxgas 2 1 1 is assumed to be the X velocity at the interface between grid cell 2 and 3 Figure XX Accordingly the first entry in the gasvelocityfile vx file will be the X velocity defined at the interface between grid cell 0 1 1 and 1 1 1 So in the case of where ContinuousRead is specified as the file format the code would read READ 52 qxgas jx jy 32 jx 0 nx jy 1 ny jz 1 nz READ 52 qygas jx jy 22 jX 1 nx jy 0 ny jz 1 nz READ 52 qzgas jx jy 32 jx 1 nx jy 1 ny jz 0 nz Calculate flow Keyword followed by true or false which indicates whether flow wi
97. le 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 A A etc contain contributions from both the heterogeneous reaction term RU and from the total 11 soluble concentrations 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 50 which arise in the case of the global implicit method should take into account both the sparseness of the Jacobian matrix A 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 These parameters are located in data statements at the top of the named routines The tolerance on function residuals is co
98. liquid partitioning gases must be explicitly entered In the case of unsaturated transport and reaction an actual mass balance on the individual gas concentrations is carried out If the database sweep option is specified then the code will automatically sweep the thermodynamic database looking for relevant gases based on the user s choice of primary species Any entries in the GASES keyword block will be ignored in this case These will then be listed in the standard output file input filename prefix out Database Format The format for gases in the database is similar to that of the secondary species lt Gas gt lt Molar volume not currently used by code Number of species in reaction integer gt Stoichiometric coefficient lt SpeciesName gt lt Stoichiometric coefficient lt SpeciesName gt lt Log K array Molecular weight Input File Entry of Secondary Species Gases are input as a simple list with one gas per line of the GASES keyword block If the database sweep option has been chosen previously the list of gases written to the file input filename prefix out may be copied from this file and then pasted directly into the GASES block MINERALS The keyword block MINERALS refers to any solid aqueous phase heterogeneous reaction mineral water reactions being the most common of these The keyword parameters for this block are simply the names of the solid phases involved in the reaction Crun
99. ll be calculated within the code or not Syntax calculate flow logical Default False Explanation This keyword indicates whether flow will be calculated according to Darcy s Law or not At this time only fully saturated flow calculations are possible Once this logical is set to true the code will look for inputs on the pressure and permeability fields Read PermeabilityFile Keyword followed by a file name containing permeability values over the entire spatial domain Syntax read PermeabilityFile filename format Backwards compatible read permeability Default None 80 CrunchFlow Manual Carl I Steefel Explanation This keyword provides a filename prefix for ASCII files containing permeability values over the entire spatial domain The extensions x y and z are assumed to indicate the files containing the X Y and Z permeabilities respectively For example the syntax read PermeabilityFile copperleach SingleColumn would open and attempt to read the files copperleach x copperleach y copperleach z using the SingleColumn see Section 7 1 2 format Units of the permeability are m The format of the file depends on the file format specified In none is provided a single column format consisting of a single permeability per line with NX varying first then NY and then NZ is assumed Unlike the fluid velocity permeabilities are defined at the center of the nodes since at least in some cases they ar
100. lls of 2 meter spacing and 10 cells of 4 meter spacing xzones 10 1 0 10 2 0 10 4 0 resulting in a total of 30 cells in the X direction Lines may be continued on the following line by using an ampersand amp at the end of the line as a continuation marker Yzones Discretization in the first Y direction Syntax yzones Zcells spacing Zcells spacing amp cells spacing where cells is an integer and spacing is a real number Default 1 0 if yzones is not specified Zzones Discretization in the first Z direction Syntax zzones cells spacing fcells spacing amp cells spacing where Zcells 1s an integer and spacing is a real number Default 7 1 0 if zzones is not specified PRIMARY SPECIES The PRIMARY SPECIES keyword block is used to specify the primary component species for a particular problem Borrowing from the tradition of EQ3 EQ6 there is a list of database primary species given in the first section of any database file The list of database primary species serves as the building blocks of reactions in the main part of the database From the point of view of CrunchFlow however primary species are any that are so designated in the input file The user may list any mathematically valid set of primary species which spans the concentration space for a particular problem If aqueous kinetic reactions are included in the problem then these must be made up entirely of the user s choice of primary species n
101. lowed by either the label zone or default If zone is specified then the code expects the X Y and Z coordinates in the form of a beginning JX and ending JX a beginning JY and an ending JY and a beginning JZ and an ending JZ in each case separated by a hyphen If default is specified instead of zone the code will initialize the entire spatial domain 1 e all grid cells to the value More than one specification of permeability_x can be and normally is provided For example in a 2D system with 20 grid cells in the X coordinate direction and 20 in the Y coordinate direction a higher X permeability zone running down the center of the domain could be specified with permeability x 1 E 13 default permeability x 1 E 12 one 0 2 10 10 1 1 Note that values of the X permeability at grid cell 0 and NX 1 have been specified since the harmonic mean of the permeability is calculated at the system boundaries using these ghost cell values Zones specified later in the sequence overwrite zones specified earlier in the example above the default specification sets the permeability x value to 1 E 13 m initially in all of the grid cells while the next permeability x keyword creates a higher permeability zone in grid cells JY 10 over the entire X length of the domain To create a no flow boundary at JX 0 over the entire Y length one could use permeability x 1 E 13 permeability x 1 E 12
102. lution max Maximum decrease in mineral volume fraction for any one time step Syntax dissolution max value value is a real positive number Default 0 001 m mineral m porous medium Explanation This keyword parameter specifies the maximum decrease in mineral volume fraction dimensionless for any one time step This option is applied in calculating the size of the time step for the next time step based on rates computed in the preceding one It is not therefore used currently as a firm control on time step since no repeat of a time step is carried out if dissolution max is exceeded in a step Fix saturation A constant uniform liquid saturation may be specified with this keyword Syntax fix saturation value value is a positive real number between 0 0 and 1 0 Default 1 0 Explanation This is an optional parameter that can be used to set a constant uniform liquid saturation for a problem The value must be greater than 0 0 and less than or equal to 1 0 Generic rates Keyword followed by a value representing the logarithm of the reaction rate to be used for all minerals in the MINERAL keyword block This option avoids the need for kinetic entries for the minerals in the database Syntax generic rates value value represents the logarithm of the mineral rate in units of mol m sec Default The option is not selected and individual kinetic entries will be required in the database Explanation The keyword parameter generic r
103. lux 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 particular 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 which is needed to run GIMRT 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 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 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 yr Diffusion coefficients are in units of m s All other references to time time for 25 CrunchFlow Manual C
104. mestep init time tolerance These keyword parameters are described in detail below Example RUNTIME time units CrunchFlow Manual Carl I Steefel timestep init 1 E 09 timestep max 0 01 time tolerance 0 05 gimrt true debye huckel true database sweep false speciate only false master H screen output 10 END Coordinate Keyword followed by an option either rectangular cylindrical or spherical Syntax coordinate option Default Rectangular Explanation Rectangular coordinates are the default but the user may also choose cylindrical or spherical coordinates Correction_max Absolute magnitude of the maximum correction to the natural logarithm of a component species in any one Newton step Syntax correction_max value value is a positive real number Default 2 0 Explanation This is an optional parameter which controls the absolute magnitude of the maximum allowable correction to a component species in any particular Newton step The correction is in units of the natural logarithm of the species molality The limitation of the correction is a form of damping of the Newton step which prevents overflows in many cases Correction factors over 10 are not recommended Excessively small correction factors lt 1 will slow the rate of convergence of the Newton method significantly Courant_number Maximum Courant or Courant Friedrichs Lewy number defined by CFL 7A TA for OS3D option GI
105. mputed from tolmax atol rtol C 58 J where C refers to the primary species concentrations If they are changed the code should be recompiled The functions residuals when multiplied by the timestep de1t in order to compare to the tolerances atol and rtol have units of moles per m bulk volume The default setting for the absolute tolerance is 10 mole m which corresponds to 10 mole liter Convergence will also be obtained where either the absolute value of the log concentration 21 CrunchFlow Manual Carl I Steefel correction or the linear concentration correction drops below 10 i e approaching the roundoff error of the computer Updating Mineral Properties When convergence is achieved the mineral volumes and surface areas and porosity if desired are updated The individual mineral volume fractions 4 are updated with the expression Q t At 7 Q t V r t At At 59 where V is molar volume of the mineral and r is its reaction rate There are now threed different options for treating porosity described fully in the section on the input file One option is to specify a single value for the porosity which then remains constant for the duration of the simulation This is set at the top of the input file using the parameters jpor and constantpor In this case the value of constantpor overrides any choice of the mineral volume fractions below The second option is to have the porosity computed from
106. ms Water Resources Res 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 1857 1873 Gupta A D Lake L W Pope G A Sephernoori K and King M J 1991 High resolution monotonic schemes for reservoir fluid flow simulation Jn Situ 289 317 Helgeson H C 1968 Evaluation of irreversible reactions in geochemical processes involving minerals and aqueous solutions I Thermodynamic relations Geochim Cosmochim Acta 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 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 1089 1198 Helgeson H C Garrels R M and MacKenzie F T 1969 Evaluation of irreversible reactions in geochemical processes involving minerals and aqueous solutions II Applications Geochim Cosmochim Acta 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 N
107. n and only periodic information is needed This number should be a minimum of 1 output every time step 41 CrunchFlow Manual Carl I Steefel Solver Keyword followed by an option either bicg gmres bcgs or direct Syntax solver option Default bicg Explanation This sets the solver method for PETSc to solve the flow or diffusion equations It does not apply to the solution of the global implicit reactive transport equations in the GIMRT option The choices include a bi conjugate gradient solver bicg and the additional options of gmres bcgs and direct The direct option carries out a direct solution of the matrix without iteration 1 e it is not a sparse matrix solver See the PETSc web pages at www anl gov petsc for further explanation Speciate only Keyword followed by true or false which selects whether or not to speciate the geochemical conditions and then stop Syntax speciate only ogical logical 1s a standard Fortran logical true or false Default false Explanation The keyword parameter speciate only is followed by true or false or yes or no and is used to select when true the option to speciate the geochemical conditions and then to stop without carrying out any reactive transport calculations This option is contrasted with database sweep which in addition to stopping after the speciation calculations also sweeps the database to load all of the possible species and minerals Using speciate only the user
108. necessary to define both a mobile total concentration 9 CrunchFlow Manual Carl I Steefel N bil bil bil Qe E Ee S ee 14 J a al and an immobile total concentration made up of the surface complexes N jenes Gime Y y xe 15 J J 1 1 The governing differential equation then becomes 0 mobile immobile ze M QU r gr mH 16 V e up M y U DV p M U7 R 9 Lave N The concept of a total concentration will be clarified with a few examples Consider a system in which the aqueous species consist of H O H OH Al and Al OH Assuming that these species are all in equilibrium then the number of unknown concentrations i e N is reduced from the total number of species concentrations by the number of independent equilibrium reactions N between them AI 4H O Al OH AH 17 H OH H 0 18 Choosing H O H and Al as the primary species and writing the chemical reactions in matrix form as in Equation 8 that is with the production or destruction of one mole of the secondary species the reactions become 4 4 4 Bop Al R teen e 1 1 0 OH which implies that the total concentrations are Uno Cno 4C rom Con O AK OH Uir C UM HORT Con 21 U C C 22 Al Al Al OH 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 speci
109. ngle quotes Input File Entry of Primary Species The format of the PRIMARY SPECIES keyword block is simple consisting of a list of species name to be used In contrast to other keyword parameters in the input file the species names are case sensitive This is to remove the ambiguity which may arise from various species and compounds have similar names for example carbon monoxide or CO and cobalt or Co The list of species in the PRIMARY SPECIES keyword block is important in that it defines the chemical space for a particular problem thus the term component species When geochemical conditions are given all of the primary species listed in the PRIMARY SPECIES keyword block must be specified in each of the conditions In this respect CrunchFlow differs from a code like PHREEQC Parkhurst 1995 where solutions may be mixed which have different sets of non zero chemical concentrations specified PHREEQC presumably sets unmentioned concentrations to 0 Since CrunchFlow grew out of a standard finite difference approach to solving partial differential equations rather than out of reaction path approaches it requires a consistent set of primary species for each cell within the domain and for the system boundaries The two examples below represent two different basis sets of primary species which can be used to solve the same problem They are mathematically equivalent therefore given the correct choice of constraints The total concentrat
110. ni LT CrunchFlow Manual Carl I Steefel LOWE DS ACTI TT C 21 Constant MOW aiia Ros etas EORR RR sosro Ua LU TAS RETIRER CAE sasso aois QR coissa soseo Fes e de Deu Nu D2 CONSCATAE DASE OW Af 1 bins 78 CAS pup os ssssrecsssi soets sesri cesis cs cnusie eisai suns tvsuesvneees oe iaai ia LO Read Velocity EMe ete te EROR D RITU I d RODA n MR dU Read Gas V Clocit y Pale e 22 Caleulate TOW seessesooesoossesoossossseseossosssessossoossessossoossossossoossossssssossosssesoossosssesssssssssessssssse OU Read Periniea Dility Ji cepe OU Permi abilily Koieiicas le cesrasceeveratszsateeacaeacadavesieusnacsaavlscueatees Ge REVO URL eueidisuecuaseateouaswaselatesiiaeae D Pe rini ali by db err TL Permeability ATTE ETSI LATTE MX PAY MEY Mc c H OO Pr ss te ire etti OF Initialize hydr static 1 ieeiriceeee ne tape aeuo rana eoa aa eran aaa us onsec TEMPERATURE seien da OD Set temperature I Temperature pradietit 5 eerte poer t i aet DER HER EH REFER IEREE PUE RM NN resos s epe a P nee Pi Ue ees OO Read T mperatu reEle iiiui une ette sap FRcERA FIG earl Dear e OD M GRo Fra sa ER RL HEROS Sd EP EET ERR PRX RUE OD BEEK VOT OSICY de OO IVEATARTAAULTEN MONETE TT DE IU ESSE aret VIR AMBITULCE TT PP OL R ad PorosityFil ous eo i tex erre pes RE eek p ere ne eui t CreatePes
111. not a sparse matrix solver See the PETSc web pages at www anl gov petsc for further explanation Pclevel Keyword followed by an integer value giving the level of fill to be used in preconditioning the flow or diffusion matrices Syntax pclevel integer value Default 5 Explanation This sets the preconditioner level of fill for solving the flow or diffusion equation An integer value must be provided Normally levels of fill of above 5 should only be used in the case of an ill conditioned matrix as typically occurs with large differences in permeability Use of higher fill levels will slow execution in the case of reasonably well conditioned linear systems See the PETSc web pages at www anl gov petsc for further explanation Reaction path Keyword followed by true or false which selects whether or not to run in reaction path mode as opposed to batch mode Syntax reaction path ogical logical 1s a standard FORTRAN logical true or false Default false Explanation The keyword parameter reaction path is followed by true or false or yes or no and is used to select when true the option to run in reaction path mode This mode is contrasted with batch mode and applies only in 0 zero spatial dimension problems a single grid cell is used In reaction path mode mineral volume fractions are not updated This is meant to represent loosely the case of pure advective transport in which secondary mineral product
112. not currently be used as a primary species in GIMRT since gases are not included in the mass balances One could use O however or any number of other species which form is used as the basis species in redox reactions in the EQ3 EQ6 data base it must be 2 aq a redox couple for example Fe and Fe or HS and SO Acid Base Reactions Acid base reactions require no special treatment other than in the initialization process where one normally fixes or calculates pH rather than an H concentration Once a pH has been either fixed or calculated in the initialization procedure however H can be treated like any other component Again no special considerations are needed in contrast to the statements by Yeh and Tripathi 1989 Activity Coefficient Model At this time the user may choose to run the code either with unit activity coefficients or with an extended Debye Hiickel formulation to calculate activity coefficients for the ionic species The Debye H ckel formulation is given by AZ I 1 2 i logy w bI 63 1 B d CU The coefficients for the temperature dependence of the parameters are taken from the EQ3 data base Initialization Procedure GIMRT initializes the solute concentrations at individual grid points over the domain 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
113. o 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 advective 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 cases 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 condition 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 example 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 O or CO is present in such large abundance 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
114. oH SEN NE EPA SEN ORNA 9E Caes ete seis eori esami a L9 Solving the Equations with Newton s Method e c eee esee ee eene eren ee eese sees seeees sees 19 Structure of the Jacobian Matrix siccscsccssccssensssessesssvennsveacsoosseoscsencsseassoassossneesestaassonsssonsesssscees 20 Convergence Criteria c esoceeseorvtet apos ox Eee CPRR Y Se Y e Rie eor eio ose re te ridi tiet e nie AL Updating Mineral Prop ertles o nee oe nean vo arp eo eoo eo ota o iato ror tu paene ren ee nene n neon ra eno 22 addu 2 ues AD Boundary Conditions use eeseeceos D cans eerie NN eIR ccoocasvaceaxeok SEP Ur VEU UR NY croce Eee E eme vreRe NV pee vas 20 Choice of Primary and Secondary Species eere eee eee e eerte ee eere ee eese eese sese es seee 27 Redox ROactlolls eerie ioci ei don teda oso aera OPER p en adn en D RN Kobe p RE Ree PAVO RM Ce DA RQe MRRRaK EK essa ve Dodd LO Acid Base Reactions ees cetis base pte bae Vlad p er PA ee HERR okoo sisses pises sosete siss iois osasse ERES QUERER ue AO Activity Coefficient Model 21 Initialization Procedure 2 Temperature Gr3dients eise eet decies ee ker ee PR e eo eee AN PRI Fei pF edv dO na neve H Pi dad PRIV da e S E dek D Rung Crunch lO E PL Reading Input File Name from a File 4 eese eee eese ee eee eee ee eee eese eo sese es seeees sees OO Keyword BIOCKS dons le sbat
115. ode i e f 0 56 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 34 The Jacobian matrices exhibit a block tridiagonal structure in a one dimensional system The blocks making up the global matrix A are N by N submatrices corresponding to the primary species at any grid point That is if we are writing the conservation equation for SiO at a particular point in space the 2 aq 20 CrunchFlow Manual Carl I Steefel AL 2 aq equation fso Will be differentiated with respect to all N primary species e g SiO 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 A A 0 0 A A A A 0 57 0 A A A 0 0 A A 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 sing
116. omplexes Ns Nexs Nsurf T E C ae VC E 2s VC m 254 vC 2 where Ck refers to the Nex exchange species and CIl refers to the Nsurf surface complexes To specify a total concentration within a condition one uses the form lt PrimarySpeciesName gt lt Concentration gt or equivalently lt PrimarySpeciesName gt lt Concentration gt total where the total is assumed to be the default To distribute the specified mass over the exchangers and surface hydroxyl sites as well the user follows the concentration value with the keyword equilibrate_surface as in the first three entries of the following example Am 1 00e 10 equilibrate_surface Np H 1 00e 10 equilibrate_surface 66 CrunchFlow Manual Carl I Steefel Put 241 1 00e 10 equilibrate surface K 8 70e 5 Mg 2 06e 5 In this example the Am Np and Pu will be partitioned between the aqueous phase and the exchange and surface hydroxyl sites 1f present while the total concentration in the case of K and Mg will be distributed completely within the aqueous phase Species Concentration Constraint The species concentration constraint is straightforward and need not be discussed here other than to show how it is specified in the input file Since the total concentration constraint is considered the default it is necessary to specify the keyword species after the concentration value lt PrimarySpeciesName gt Concentration species
117. on expressed as log Q Keq where the dissolution rates 1 2 of its maximum value However normally it is preferred to use the DissolutionOnly in a kinetic entry in the database For the case of linear kinetics the code uses an analytical expression for the approach to equilibrium that is continuously differentiable with the rate given by E R 20 5 0 75 4 0 25 p its E 57 where the quantity is defined by CrunchFlow Manual Carl I Steefel For nonlinear kinetics the factor of is not used and precipitation is simply suppressed while the rate law provided in the database is used to give the dependence on Gibbs free energy Local equilibrium with respect to mineral phases variously known as the local equilibrium approximation or fantasy depending on one s point of view can be specified through a combination of parameters To suppress the reduction in reactive surface area as the mineral volume fraction of gt 0 the following option may be specified Reaction name local equilibrium In other words the initial reactive surface area remains fixed even during dissolution as long as the mineral is present at all This is of course unrealistic but the local equilibrium fantasy has been so drummed into a generation of geochemists and petrologists that it is impossible to dispense with this option This alone does not enforce local equilibrium but when combined with a sufficiently large rate constant and reactive s
118. one or more selectively with something like default 0 0 1 20 1 1 permeability x 1 E 13 default permeability x 1 E 12 on 0 0 1 8 1 permeability x 1 E 12 on 0 0 12 20 1 1 which would then create a no flow boundary at JX 0 except for JY nodes 9 10 and 11 The convention adopted here is that all three dimensions must be specified as zones even in a one dimensional problem although permeability y and permeability z are not needed 82 CrunchFlow Manual Carl I Steefel Permeability_y Keyword used to specify the permeability in the Y coordinate direction by zone in the input file Syntax permeability y value zone jxbegin jxend jybegin jyend jzbegin jzend or permeability_y value default Default None Explanation The permeability_y keyword uses an input format similar to other parameters specified by zone e g pressure tortuosity See the keyword permeability_x for details on how this keyword is used to specify permeabilities by zone The convention adopted here is that all three dimensions must be specified even in a one dimensional problem Permeability_z Keyword used to specify the permeability in the Z coordinate direction by zone in the input file Syntax permeability_z value zone jxbegin jxend jybegin jyend jzbegin jzend or permeability_z value default Default None Explanation The permeability_z keyword uses an input format similar to other p
119. ons of the individual species at the end of the first 1 2 transport step U are the total concentrations at the end of the full reaction step and C7 are the individual species concentrations at the end of the second 1 2 transport step 1 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 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 ww lt 4 37 Ot 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 error The global implicit approach of course eliminates this error but at the cost of adding in transport errors associated with the spatial discretization 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 approach Our experience however is that the sequential iteration scheme offers little improvement over the Strang operating splitting procedure and that
120. ot necessarily those found in the database It should be noted however that some choices of primary species will result in non intuitive components making it difficult to specify initial and boundary conditions for a particular problem The primary component species or basis species are the chemical building blocks for a particular problem The concentrations of the primary species are the major unknowns in the chemistry part of a reactive transport problem The partial differential equations for the conservation of chemical mass in the system are written in terms of the total concentration of these primary component species This formulation which has been widely discussed by a number of workers is used instead of equations for the conservation of the various chemical elements 47 CrunchFlow Manual Carl I Steefel Database Format Database primary species are listed in the database immediately following the specification of the temperature field first line and 3 lines giving the temperature dependent Debye Huckel parameters The primary species block in the database is terminated with the line End of primary 0 0 0 0 0 0 The database primary species are specified in the following form lt SpeciesName gt lt Debye Huckel size parameter gt lt Charge gt lt Molecular weight gt Each database primary species is given on a separate line in the database in free format The species name is enclosed in si
121. ous kinetics Entries in the aqueous kinetics section of the database take the form lt Aqueous reaction label gt Number of parallel rate laws lt ReactionType gt Number of species in reaction integer lt Stoichiometric coefficient lt SpeciesName gt lt Stoichiometric coefficient lt SpeciesName gt Log K gt which is then followed by the number of lines corresponding to the number of parallel rate laws which are described below A minimum of one rate law needs to be given in the database following the kinetic reaction stoichiometry Reaction type specified in the first line refers to rate laws of the type tst monod or irreversible just as in the solid liquid phase reaction input although the format of the input clearly differs A fuller description of the various reaction types is given below Types of Rate Laws TST The tst rate law is described in detail by Lasaga 1981 1984 and by Aagaard and Helgeson 1981 and takes the form where ks is the intrinsic rate constant in units of mol kg water yr Q is the ion activity product for the mineral water reaction Keg is the corresponding equilibrium constant and le is a product representing the inhibition or catalysis of the reaction by various ions in solution raised to the power n Unlike the database entries for solid liquid phase reactions there is no dependence on surface area and currently no allowance for a temperature dep
122. ous medium In this way a mineral which is neglected in a particular condition is treated as a potential secondary phase that is it is not initially present but it could form if supersaturated Since boundary conditions and source terms represent ghost cells or zones outside of the domain which are therefore not computed as a function of time supersaturation of a mineral phase will have no effect Mineral Volume Fractions The first values to be input following a mineral solid phase name is the volume fraction The ability to add mineral concentrations in other units mass fraction etc will be added in future releases The format is name of solid phase volume fraction Surface area information is appended to the same line and is described below Mineral Surface Area Options Mineral surface area is important in two ways 1 it gives the reactive surface area to be used in mineral dissolution and precipitation reactions and 2 it provides the surface area upon which adsorption surface complexation may occur see section on Surface Complexation Input There are currently two options for specifying mineral surface area 1 bulk surface area in units of m solid phase m porous medium or 2 specific surface area in units of m g The format is name of solid phase volume fraction surface area option value where the surface area options are bulk surface area or bsa and specific surface area or ssa respectiv
123. ous species in solution dependence lt species name gt lt exponent gt species name lt exponent gt The dependence of the far from equilibrium rate is given by pairs of species names in free format followed by the appropriate exponent An example of a database input which provides both pH dependent and pH independent calcite dissolution rates is given below UELLE eee Oe NE Tee AE Ae ce em Calcite label default type tst rate 25C 6 19 activation 15 0 kcal mole dependence Calcite label h 53 CrunchFlow Manual Carl I Steefel type tst rate 25C 3 00 activation 15 0 kcal mole dependence H 1 0 A general form for a dependence on reaction affinity or Gibbs free energy is given following that given in Burch et al 1993 and Hellmann and Tisserand 2006 The affinity dependence of the rate is defined with the parameters m m and m3 A more familiar form would be rn d Ka AffinityDependence ml m2 m3 The form in the database is As an example the rate law for albite dissolution proposed by Hellmann and Tisserand 2006 can be implemented by specifying two parallel rate laws Albite label HellmannFFE type tst rate 25C 11 9897 activation 0 00 kcal mol dependence AffinityDependence 1 00 0 0000798 3 81 Albite label HellmannCTE type tst rate 25C 13 743 activation 0 00 kcal mol dependence AffinityDependence 1 17 1 00
124. parameters may appear anywhere within a particular block but certain keyword parameters must appear in the appropriate blocks The possible keyword blocks include TITLE PEST RUNTIME DISCRETIZATION OUTPUT INITIAL_CONDITIONS PRIMARY_SPECIES BOUNDARY_CONDITIONS SECONDARY_SPECIES TRANSPORT GASES FLOW MINERALS POROSITY AQUEOUS_KINETICS TEMPERATURE ION_EXCHANGE CONDITION SURFACE_COMPLEXATION With the exception of the keyword block CONDITION the blocks should appear only once in the input file Each keyword block is terminated by an END The keyword block CONDITION is a special case in that it can occur multiple times Each occurrence of CONDITION specifies a separate geochemical condition these may be boundary or initial conditions or source terms containing the geochemical input needed to describe a particular problem The various keyword blocks will be discussed individually in more detail below Reading Input File Name from a File Normally CrunchFlow will prompt the user to enter the name of an input file interactively The user can provide the name of the input file however by including it in a file names PestControl ant The code will check to see if this file exists in the directly from which CrunchFlow is invoked and if it does will attempt to read the file name from it Upon successfully reading the input file name the code will then check to see that this file exists before reading from it The use of the PestControl ant to p
125. possible spatial units which can be specified with the space units keyword parameter Time Units Alternate Acceptable Specifications years default year days day hours hour minutes minute seconds second Table 2 List of possible time units which can be specified with the time units keyword parameter 31 CrunchFlow Manual Carl I Steefel Input File Formats Various spatially dependent field variables can be calculated outside of CrunchFlow and read in for use These include the temperature permeability tortuosity porosity liquid saturation and erosion burial rates Typically these are specified with a keyword see below followed by the file name As of August 2006 it is also possible to specify the format of the file The basic syntax is keyword filename InputFileFormat where the keyword may be one of the following read VelocityFile read GasVelocityFile read PermeabilityFile read PorosityFile read SaturationFile read TortuosityFile read BurialFile read TemperatureFile The possible file formats include SingleColumn the default ContinuousRead FullFormat dummy X and optionally Y and Z coordinates Unformatted binary SingleColumn Example DOjz 1 nz DO jy Lny DO jx 1 nx READ 52 tortuosity jx jy jz END DO END DO END DO ContinuousRead Example EAD 52 qx jx jy jz jx 70 nx jy 1 ny jz 1 nz EAD 53 qy jx Jy J32Z jJx 1 nx jJy 0 ny jz 1 nz EAD 54
126. quivalently 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 asecond order tensor which is defined as the sum of the mechanical or kinematic dispersion coefficient D and the molecular diffusion coefficient D in water divided by the formation resistivity factor F D D D 4 ze 0 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 example Dullien 1979 but here a definition based on Archie s Law is used which gives the formation factor as F po 5 where m is the cementation exponent Dullien 1979 reports values of m which range between about 1 3 and 2 5 The kinematic dispersion coefficient is written here as 2 Uu D a u a m e 6 u CrunchFlow Manual Carl I Steefel where a and a are the longitudinal and transverse dispersivities respectively Bear 1979 Marsily 1986 and u is the magnitude of the velocity Off diagonal dispersion coefficients i e D j Where i and j are coordinate directions are not included in GIMRT at the present time The heterogeneous reaction term R can be written as the sum of all the individual mineral water reactions which affect the concentration of the 2 th species Nin R
127. r advection resulting from variations in the velocity field that cannot be resolved by the numerical grid In this standard representation based on Fick s law dispersion is driven by the gradient of concentration and parameterized by a linear proportionality constant 1 e the diffusion coefficient This is a gross simplification of a complex and poorly understood phenomenon As such the formulation neglects many features 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 material 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 Solving the Equations with Newton s Method As an example consider a problem involving one dimensional multicomponent reactive transport Equation 49 for any component at the grid point P is a nonlinear function of the concentrations of the primary species at the P E and W nodal points using the global implicit approach of GIMRT An iterative scheme therefore is required to solve the set of nonlinear equations in either case The code solves the nonlinear equations
128. racked as a time series The time series print keyword may be followed with all in which case all of the total concentrations will be written out in column form in units of mol kg followed by all of the species concentrations in logarithmic form Listing of species names will cause the total concentrations of these species to be written out not the individual species concentrations In addition pH may be specified to be output As an example the following line time series print pH Cat Na K will cause the code to track the time evolution of the solution pH and the total concentrations of Ca Na and K Species names must match those in the primary species list case sensitive 45 CrunchFlow Manual Carl I Steefel Time_series_interval Time step interval at which a time series should be output Syntax time series interval interval where interval is an integer Default 1 Explanation This keyword parameter gives the time step interval at which the time evolution of the species specified in time series print will be written out at the nodal location specified by time series keyword Time series output Specific times at which a time series should be output Syntax time series output time time2 time3 where timel and time2 are real numbers Default None Explanation This keyword parameter gives the option of writing out time series at specific times This is useful when the objective is to compare discrete dat
129. racters contining the liquid saturation distributed over the entire spatial domain Default None Explanation This provides the name of a file containing temperatures distributed in space to be read This option supersedes all other temperature specifications The format of the file depends on the optional format specified see Section 7 1 2 If no format is specified the code assumes a single column SingleColumn of values in the order 1 NX 1 NY 1 NZ DO jz 1 nz DO jy 1 ny DO jx 1 nx READ 52 t jx jy jz END DO END DO END DO POROSITY This keyword block is used to set parameters and options related to the treatment of the porosity The main contrast is between simulations in which the porosity is fixed and those in which it evolves as the mineral volume fractions evolve Fix_porosity Keyword followed by a value giving the default fixed porosity Syntax fix_porosity value Default If not provided the code calculates porosity based on the mineral volume fractions Explanation The keyword parameter fix_porosity is used to set a global fixed porosity for a simulation When provided it disables the calculation of porosity based on the sum of the volume fractions Accordingly no porosity update based on these quantities is possible When this parameter is set in the POROSITY keyword block set_porosity keywords within individual geochemical conditions can be used once the geochemical condition is distribu
130. re defined at grid interfaces not at the center of the nodes so NX 1 X velocities are needed in the X coordinate direction NY 1 Y velocities in the Y coordinate direction and NZ 1 Z velocities in the Z coordinate direction CrunchFlow uses the convention that velocities for a particular grid cell are defined at the right hand interface the so called right hand rule So for example the velocity qx 2 1 1 is assumed to be the X velocity at the interface between grid cell 2 and 3 Figure XX Accordingly the first entry in the velocityfile vx file will be the X velocity defined at the interface between grid cell 0 1 1 and 1 1 1 So in the case of where ContinuousRead is specified as the file format the code would read READ 52 qx jx jY 3Z jx 0 nx jy 1 ny jz 1 nz READ 53 qy jx 3Y 3Z jx 1 nx jy 0 ny jz 1 nz READ 54 qz jx jY 3Z jX 1 nx jy 1 ny jz 0 nz Read GasVelocityFile Keyword followed by a name for a file containing gas velocities over the entire spatial domain Syntax read GasVelocityFile filename format Backwards compatible read gasvelocity Default None Explanation This keyword provides a filename prefix for ASCII files containing gas velocities over the entire spatial domain The extensions vx vy and vz are assumed to 79 CrunchFlow Manual Carl I Steefel indicate the files containing the X Y and Z velocities respectively For example the syntax read GasVelocityFi
131. rn off the database sweep option The list of secondary species generated by the database sweep are found in the output file input filename prefix out which is generated every time a run with CrunchFlow is carried out This is done because of the large number of potentially extraneous species which may be loaded for a particular problem especially in the case of organics The disadvantage of the approach is that the user may potentially neglect some important secondary species 1 e those which have non trivial concentrations in a particular problem The recommended procedure is to begin a problem by running the code in database sweep mode and then winnowing the list of secondary species if necessary In most one dimensional problems the full list of secondary species may be used with the exception of problems where a huge list of organic species may be loaded The format for listing secondary species in the input file is the same as the format for primary species Secondary species to be included are listed on separate lines within the SECONDARY SPECIES keyword block If the database sweep option has been chosen previously secondary species written to the file input filename prefix out may be copied and then pasted into the input file As with primary species the input is case sensitive the species listed must be found in the database or an error message will result Example ECONDARY SPECIES GASES T
132. rovide an input file name will then 30 CrunchFlow Manual Carl I Steefel skip the requirement of interactive input from the user an option that is particularly useful when running PEST Keyword Blocks Each keyword block is initiated with the appropriate keyword case insensitive on a line of its own and terminated by END also on a line of its own The contents of any keyword block then is anything between the keyword block name and the END specifier Blank lines within the keyword block will be ignored Lines beginning with are treated as comments and ignored Space and Time Units Keyword parameters involving space and time units are used in several of the keyword blocks These take the form time_units unit space units unit where unit in the case of time may be years the default days hours minutes or seconds and in the case of space may be meters the default kilometers centimeters millimeters and micrometers Non metric units e g feet are not recognized As with any other keyword parameter they may occur anywhere within a particular block to which they apply A specification of a space or time unit however is not global and applies only to the keyword block within which it occurs Space Units Alternate Acceptable Specifications meters default meter m kilometers kilometer km centimeters centimeter cm millimeters millimeter mm micrometers micrometer micron microns um Table 1 List of
133. s are left behind in a flow path Read SaturationFile Option to specify a file with liquid saturations to be read Syntax read SaturationFile filename format filename gives the name of the file up to 132 characters contining the liquid saturation distributed over the entire spatial domain Default None Explanation This provides the name of a file containing liquid saturations distributed in space to be read This option overrides the fix saturation keyword The format of the file depends on the optional format specified see Section 7 1 2 If no format is specified the code assumes a single column SingleColumn of values in the order 1 NX 1 NY 1 NZ DO jz 1 nz DO jy 1 ny DO jx 1 nx 40 CrunchFlow Manual Carl I Steefel READ 52 satliq jx jy jz END DO END DO END DO Restart Option to use a restart file Syntax restart filename append filename is an optional string up to 132 characters indicating the name of the restart file Default crunch rst Explanation This keyword parameter specifies that the code should be restarted If restart is not followed by anything the default restart file name of crunch rst is used Alternatively the user may specify an optional restart file name this can be useful since the default file name will be overwritten with a new run This optional restart file however must have been created in an earlier run with the save restart keyword Binary restart
134. s in 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 Figure 2 Program structure for GIMRT which uses a global implicit or one step method used to solve the multi component reaction transport equation 24 CrunchFlow Manual Carl I Steefel Specify initial and boundary conditions Begin time stepping Advect all species using TVD algorithm Diffuse and or disperse species Sequential Iteration Add source term Reaction at individual space points Calculate secondary species concentrations Calculate total concentrations Calculate mineral reaction rates Calculate Jacobian matrix Solve system of equations Update concentrations until converged Sequential Iteration 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 Features of CrunchFlow Before running the code the user of GIMRT needs to decide what physical process es he or she wants to model and design the simulation accordingly This involves choosing the boundary conditions for the problem e g no f
135. s list of secondary species and minerals given in the input file is used Steady state Keyword followed by true or false which selects whether or not to run to a quasi steady state as a condition for stopping a particular simulation Syntax steady state logical tolerance logical 1s a standard FORTRAN logical true or false tolerance is an optional real number giving the criteria for attainment of a steady state Default false with a tolerance of 1 E 07 in units of mol kg year normalized to the primary species total concentration so tolerance has units of 1 year Explanation The keyword parameter steady state is followed by true or false or yes or no and is used to select when true the option to run to a quasi steady state and then terminate the program rather than to rely on the specific output times given in the OUTPUT keyword block If true the logical specifier can be followed by an optional convergence criteria in units of normalized mol kg year for all of the primary species total concentrations the normalization then results in a convergence criteria with units of l year Depending on the value of the convergence criteria a value of 1 E 07 is the default a quasi steady state may be attained since in many cases slow dissolution of 42 CrunchFlow Manual Carl I Steefel primary mineral phases may prevent a true steady state from being achieved In many cases however it is possible to attain a quasi stationary sta
136. s per m porous medium specific C P os Aoi MW sites V The format for inputting the molar density of sites is mol m Surface complex Value Exchanger Concentrations The concentration of ions on exchangers can be set in several ways with the most important distinction being between 1 exchange on the bulk material and 2 exchange on specific minerals In the case of exchange on a specific mineral the code will link the concentration of exchange sites 1 e the cation exchange capacity or CEC to the concentration of that mineral Thus an increase in the volume fraction of a secondary phase e g clay as a result of chemical weathering or Fe oxyhydroxide as a result of corrosion will increase the cation exchange capacity of the bulk material Exchange on bulk material This option is used where no specific mineral has been identified as the exchanger phase using the on mineral syntax in the ION EXCHANGE keyword block To calculate total concentration of exchange sites per bulk porous medium the cation exchange capacity represented by that site needs to be specified along with the density and the total volume fraction of the bulk solid phase This information must be provided in each geochemical To specify the cation exchange capacity associated with each exchanger in a geochemical condition the following syntax is used Exchanger cec Value where the CEC is given in units of equivalents per gram solid T
137. speciate only or database sweep keywords set to true in the RUNTIME keyword block is it possible to proceed without specifying initial conditions for the entire spatial domain Required keyword block followed by a specification of the name of the geochemical condition to be used followed in turn by the grid cells where the initial condition is to be set The initial condition can be fixed for the entire course of the simulation using the optional appended keyword fix Syntax condition name JX JX JY JY JZ JZ fix Default None Explanation Using this keyword block geochemical conditions are distributed over the spatial domain Multiple specifications can occur Initial conditions listed later in the input file will overwrite conditions specified above The optional parameter fix can be appended to the output so as to fix this geochemical condition for the entire course of the simulation The code now assumes that if fix appears it will occur immediately after the specification of the X coordinates JX in a 1D problem immediately after the Y coordinates JY in a 2D problem and after the Z coordinate JZ in the case of a 3D problem In the example below the code will verify that the geochemical condition names quartz zone portlandite zone and nonreactive zone have been provided in the input file This keyword block therefore is read after the geochemical conditions specified in the input file are read The code will also verify that
138. system of linearized equations Once convergence is achieved 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 computes 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 23 CrunchFlow Manual Carl I Steefel 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 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 N by N 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 differ
139. t The approach since it involves moving the reaction term to the right hand side of the equations as a source term sometimes leads to temporary convergence problems but these have not been particularly severe so far In addition 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 small amounts of a mineral are present than it is possible that dissolution of the entire amount present within a grid volume may occur 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 22 CrunchFlow Manual Carl I Steefel 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 2 Fal dissolution A h Pno A 1 m m 0 2 3 6 H precipitation h where 6 is the initial volume fraction of the mineral m and 9 is the initial porosity of the medium This formulation ensures that as the volume fract
140. tInstructionFile s scscecssecsecscosessscessseveasesevessesvessunsnsesvesessaevvdonsevenstssveenssesusescesseees 7 CreatePestExchangeFile seesseese SS Exchange eeesssecesssooesssocesssocese RUNE OO EROSION essees OO Read BurialFile iso en ORIS NUBE ERA REIN KM C RUE INI ACKNOWIEU SEMEN PERRA Literature COD e ON nO ecce CrunchFlow Manual Carl I Steefel Mathematical Formulation Conservation Equations The governing equation for the conservation of solute mass is given by Steefel 1992 Steefel and Lasaga 1994 SM noC V DV o M 4C P M oC J R 2 1 2 N 1 where C is the molal concentration of a species in solution in units of moles per kg H O M jo is the mass fraction of H O p is the fluid density u is the Darcian flux D is the combined dispersion 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 N 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 such that R R R R 2 The reaction rate is written here in terms of a unit volume of rock but could be written e
141. te see Lichtner 1988 in which aqueous concentrations are nearly unchanging while mineral volume fractions are slowly evolving Timestep_max Indicates maximum time step allowed in problem Syntax timestep_max value value is a positive real number Default 1 0 year Explanation The default units are years with the option to reset the time units to one of those discussed under Time Units e g seconds days etc using the time_units parameter keyword This is the maximum allowed time step but a smaller maximum may be used if dictated by Courant number restrictions as in the case of CrunchFlow run in OS3D mode time splitting of transport and reaction Timestep_init Indicates initial time step to be used in simulation Syntax timestep_init value value is a positive real number Default 1 E 09 years Explanation The default units are years with the option to reset the time units to one of those discussed under Time Units e g seconds days etc using the time_units parameter keyword The initial time step also serves as the minimum time step which is only superseded by a restriction on the Courant Friedrich Lewy CFL number Time_tolerance Tolerance for the second derivative of the change in the natural logarithm of the master species with respect to time Syntax time_tolerance value value is a positive real number Default 0 001 Explanation Tolerance applies to the second derivative of the natural logarithm of th
142. ted in space within the INITIAL_CONDITION keyword block Minimum_porosity Keyword giving the minimum porosity value that is allowed as the mineral volume fractions evolve Syntax minimum_porosity value Default 107 Explanation This keyword is only used to set a minimum to which the porosity can evolve as the mineral volume fractions change This option is only used where porosity update is true and neither fix porosity nor read porosity are set 86 CrunchFlow Manual Carl I Steefel Porosity_update Keyword followed by a logical true or false used to determine whether the porosity is updated as the mineral volume fractions evolve Syntax porosity_update logical Default False Explanation This keyword is only used where neither the fix_porosity or read_porosity keywords have been set When true it instructs the code to update the porosity as the mineral volume fractions evolve according to Y1 4 where N is the number of minerals in the system and 4 is the volume fraction of the individual mineral Read PorosityFile Keyword followed by a file name containing permeability values over the entire spatial domain Syntax read PorosityFile filename format Backwards compatible read porosity Default None Explanation This keyword provides the name of file containing porosity values defined at the center of the grid cell for the entire spatial domain The format of the file depends on the file format specified
143. the activation energy to 0 kcal mole There is also an option for directly inputting the diffusion coefficient at temperature 29 CrunchFlow Manual Carl I Steefel 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 of the temperature field match those specified in the code The positions in the code where a user should place a routine for calculating or inputting a temperature are marked in the routine gimrt f Running CrunchFlow CrunchFlow reads a user provided input file on startup which provides the necessary physical and chemical parameters needed for a run The input file the name of which is specified by the user is keyword based so that the order of appearance does not matter In many cases it is possible to leave various optional parameters out altogether and thus allowing the code to use default parameters Keywords are grouped broadly into keyword blocks which in turn include a variety of keyword parameters Keyword blocks may appear anywhere in the input file and keyword
144. the specification of the reaction the mineral kinetic options that can be specified in the input file include label lt label gt rate rate activation activation energy suppress precipitation e local equilibrium associate mineral name To load solid liquid phase reactions with labels other than default one uses the form Reaction name label label where label is a string which must match one in the kinetic database Similarly one can overwrite the reaction rate in the database with a specification in the input file Reaction name label label rate rate where the rate is assumed to be the logarithm of the rate in units of moles m7 s at 25C The activation energy may also be overwritten Reaction name label label activation activation energy where the activation energy is in units of kcal mol Previously it was possible to specify a threshold for nucleation but implemented in the simple fashion in which it was numerical convergence was poor due to the discontinuity in rates across the nucleation threshold The Newton method for nonlinear equations performs poorly with functions that are not continuously differentiable Now one can achieve at least suppression of precipitation by using the suppress precipitation option Reaction name suppress precipitation e where e is the value of the absolute value of the undersaturati
145. the sum of the mineral volume fractions which may vary from grid cell to grid cell but to then leave the porosity constant through the course of the simulation The third option is to compute the porosity locally from the sum of the mineral volume fractions and to update it as a result of mineral precipitation and dissolution Where the mineral volume fractions are used the porosity is computed from the g 1 gt 60 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 fractions of the minerals usually evolve much more slowly than the solute concentrations However in certain cases the approach may lead to inaccurate solutions In previous versions this was dealt with in two ways 1 reducing the time step where the mineral volume fraction would go below zero and 2 limiting the amount of dissolution and or precipitation that can occur in any one time step Currently the second method is still used but now the code automatically checks to see if the computed reaction rate exceeds the available mineral within the elemental volume If it does it calculates the reaction rate needed to cause the mineral to disappear and uses this as the reaction rate rather than the initially calculated one In this way there is no problem of a mass balance error i e dissolving more mineral than is actually presen
146. urface area equilibration within about 10 log Q Kegq units of equilibrium is possible The rate constant and or reactive surface area should not be set arbitrarily high however since this leads to poor numerical convergence due to the stiffness of the governing equations A log rate constant in the range of 3 00 to 5 00 should provide local equilibrium behavior although this depends in detail on the magnitude of the transport rates The following example should provide local equilibrium with respect to the minerals portlandite and Ca oxalate even under relatively high flow rates MINERALS Portlandite rate 3 00 local equilibrium Ca Oxalate rate 3 00 local equilibrium END Another option that has been recently added is given by Reaction name associate mineral name which is used to associate more than one reaction stoichiometry and equilibrium constant with the same mineral This is relevant normally when aqueous kinetics are involved and the reaction for the same mineral is written as a kinetic reaction involving different sets of primary species The most familiar example would be the oxidation of organic carbon by various electron acceptors none of which are in equilibrium with each other and therefore they are listed as Primary Species see below One example is given by the reductive dissolution of amorphous iron hydroxide by two different mechanisms 1 the reductive dissolution mediated by the bacteria G
147. utfiles Option to specify additional input files to be read on restart Syntax later inputfiles filename filename2 filename3 filename filename2 are strings up to 132 characters indicating the name of the input files to be read following restart Default No additional input files are read 38 CrunchFlow Manual Carl I Steefel Explanation This option allows additional input files to be read automatically restart Normally this option is used in conjunction with the keywords save_restart and restart which generate and read respectively restart of pickup files that can be used to read in the spatial distribution of species and minerals from a previous run changing a boundary condition There is no a a priori restriction on the number of files The later_inputfiles option should be specified only in the first input file in the series to be read later specifications will be ignored Master Master variable to be used for time step control Syntax master name name is a string up to 132 characters representing the name of a species Default O2 aq if present otherwise H if present otherwise species number 1 Explanation The keyword parameter time tolerance is applied to this species Normally for redox problems O2 aq is used since it is the most sensitive master variable For non redox multicomponent problems H is normally used The name of the species is case sensitive unlike other parameters sp
148. ution on the far from equilibrium dissolution rate This is most commonly the solution pH or hydroxyl ion activity but may include other electrolytes as well e g Na and Cl Dove 1995 As an example consider the overall rate law for kaolinite Nagy 1995 which includes H and OH dependent dissolution rates Baoi Aaa Kaags ait J1 Sas 29 Vg oH oH kaol Note also that Equation 28 is fully compatible with the thermodynamics of the system since if the rate constant is chosen large enough relative to the transport rates in the system local equilibrium or transport controlled behavior is obtained The temperature dependence of the reaction rate constant can be expressed reasonably well via an Arrhenius 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 E 1 1 k k a G0 se R E xis 60 where E is the activation energy k is the rate constant at 25 C R is the gas constant and T is temperature in the Kelvin 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 Depending on how mineral surface area is computed it may be necessary to include a dependence of the surface area on the matrix and fracture porosity such that A 0 as the porosity 0
149. vel are then Ax E W min a Ur An UT aU t Ar asU t At C aU t At R7 t At Ax 0 49 where At is the time step and f 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 Equation 34 and are expressed in terms of the total concentrations here in order to be concise Transport Algorithm for Operator Splitting Approach 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 transport and reactions e g GIMRT Decoupling 1 e time splitting the transport algorithm from the geochemical 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 inherent 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
150. y below the porosity threshold set in the keyword threshold_porosity Syntax tortuosity_below value Default 0 0 only used where threshold_porosity is set Explanation This keyword sets the value of the tortuosity to be used below the threshold porosity This keyword is only used if threshold porosity is set Tortuosity above Keyword followed by a value giving the tortuosity above the porosity threshold set in the keyword threshold porosity Syntax tortuosity above value Default 0 0 only used where threshold porosity is set Explanation This keyword sets the value of the tortuosity to be used above the threshold porosity This keyword is only used if threshold porosity is set FLOW Constant flow Keyword followed by a value s giving the X Y and Z components of the Darcy flux Syntax constant flow qx qy qz Default 0 0 m m yr Explanation This keyword is used to set a constant Darcy flux As a vector it includes X Y and Z components In a one dimensional problem in X the Y and Z 77 CrunchFlow Manual Carl I Steefel velocities are optional Units are changed inside the FLOW block with the space_units and time_units keywords Constant_gasflow Keyword followed by a value s giving the X Y and Z components of the gas flux Syntax constant gasflow qxgas lt qygas gt lt qzgas gt Default 0 0 m m yr Explanation This keyword is used to set a constant gas flux As a vector it includes X Y

Download Pdf Manuals

Related Search

Related Contents

Franke CRX 651  Owner`s Manual  ECC99 TUBE - ThaiHdBox  Skpad SKP-FLIP-IBO3    Battery Charger  Informationen zu NetVault Backup  MAYCO DEL SURESTE, S.A. de C.V.  

Copyright © All rights reserved. Failed to retrieve file