Home

APPENDIX - Index of

1. 0 If we make the assumption of small phreatic surface slopes Bear 1972 p 259 i e IVA lt 1 the above condition can be approximated as be 5V2 N 2 9 eff OE 2 2 Flow in the unsaturated zone For the flow of a single fluid phase in the unsaturated zone we employ the module usip of NUFT It solves the mass balance equation for a single com ponent in a single fluid phase The mass balance equation for S lt 0 is as o F 5 V K V us 2 2 10 vi APPENDIX where the primary variable used is saturation S The pressure head Pw Pwg With atmospheric pressure as datum is given as a function of satu ration by the moisture retention relationship p For 0 the mass balance equation that is solved is ga V K V 7 2 2 11 with w has the primary variable The actual implementation of the switching between primary variables is described later For the relationships p and K K the program employs the Van Genuchten relationships and respectively replacing Swr A and C by S a and m respectively 2 3 Transport of a single component in a single fluid phase The module usic of NUFT is used to solve this class of problems The flow field can be solved by the module usip or the module usnt which is described later The mass balance equation for a component in a fluid phase in multiphase flow with p p pw and w7 p7 pw is omitti
2. a e pattern ex matches all strings that begin wi e char Th tt tch ll strings that begi ith the ch acters ex b The pattern ex b2 z matches all strings that begin with the characters ex followed by any number including zero of strings which are then followed by the string b2 and which end with the character z XXV1 APPENDIX c The pattern r2 xay matches all strings that begin with r2 fol lowed by any single character and then followed by the characters xay g A set of elements in the input file can be specified by a statement of the form range lt elem range0 gt lt elem range1 gt where lt elem range0 gt lt elem range1 gt are pattern strings of the element names Element names are of the form lt zone gt lt i gt lt j gt lt k gt where lt i gt lt j gt lt k gt are the indices 1 nx 7 l ny k 1 nz of the grid where nx ny nz are the number of elements in the direction of the x y and z axes Examples are that e denotes all elements which start with the letter e and that 2 3 denotes all elements that have indices i 2 and j 3 Another way that a set of elements can be specified is through the statement C index range lt imin gt lt imax gt lt jmin gt lt jmax gt lt kmin gt lt kmax gt lt imin gt lt imax gt lt jmin gt lt jmax gt lt kmin gt lt kmax gt index range which defines a union of rectangular regions each of whi
3. or log10Kz Vmag lt phase gt Mass fraction Mole fraction Porosity Saturated permeabilities m m s for usip ucsat Log base 10 of saturated permeabilities log 10 m for usnt and log 10 m s for usip ucsat Magnitude of flow velocity vector m s lvii Vel lt phase gt The following are connection variables are defined between every pair of connected elements usually those between neighboring elements Specific fluxes are defined to be the total flux between elements divided by the flow area between the elements The units of flux are given below for the usnt module For ucsat and usip the units of flux are volumetric that is m s for specific flux synonymous with specific discharge and m s for total flux Connection Variables V lt phase gt q lt phase gt q lt comp gt q lt comp gt lt phase gt qcond Q lt phase gt Q lt comp gt Q lt comp gt lt phase gt Qcond APPENDIX Components of the flow velocity vector in the x y z directions at each element centroid valid only for field format contour or contour output options Flow velocity m s Specific flux of phase kg s m Specific flux of component kg s m7 in all phases if lt comp gt is energy then specific energy flux W m in all phases including solid phase is output Specific flux of component kg s m in a particular phase if lt comp gt is energy then specific ene
4. setcomp lt comp gt table lt t0 gt lt x0 gt lt t1 gt lt x1 gt Specify concentration wae 3 of each component lt comp gt table lt tO gt lt x0 gt lt t1i gt lt x1 gt enthalpy lt tO gt lt H1 gt lt t1 gt lt H1 gt Specify enthalpy of setcomp specified phase not present if isothermal flux allocation options see below phasef lux 3 End of the source term One can have any number of additional phase fluxes compf lux Present if one wishes to specify a 33 given flux of a component name lt source name gt 3 Name of source term used in flux output option comp lt comp name gt Name of component for which flux is 4 USING THE COMPUTER PROGRAM xlv specified range lt range gt lt range gt Range of elements table lt t0 gt lt q0 gt lt t1 gt lt qi gt Table of flux of the component enthalpy lt tO gt lt H1 gt lt t1 gt lt H1 gt Table of enthalpy of the component flux allocation options see below compf lux End of the source term One can have any number of component fluxes srctab 3 End of the all source terms Note that fluxes can be specified either as a phase flux that is a specified flux of a phase or as as component flux which is a specified flux of a given component The elements in the specified range can be but do not have to be part of the bc 2nd type command
5. 2 20 where constants A B and C are specified by the user In addition to the mass balance equations under nonisothermal condi tions NUFT solves the energy balance equation Ji el PSaPatta 1 6 pstpa T To a _V bj X bS alo rv 5 rozwa A VT a 7 2 21 where ug and ha are the internal energy and the enthalpy respectively of the a phase cps is the heat capacity of the solid matrix T is a reference temperature D is the coefficient od hydrodynamic dispersion of the y component in the a phase h are partial enthalpies and A X Aa T D a is the combined thermal conductivity and coefficient of thermal dispersion of the porous medium as a whole Since D is a function of flow velocity at this time NUFT does not have velocity dependent thermal dispersion coefficient The coefficient A e g for a liquid gas system is given by either a s VM Sel Acts s Sg Agts s 3 NUMERICAL MODEL ix or Aa s Sel Arts s Sy Agts s NUFT also has hard wired correlations for aqueous phase viscosity and densities and specific enthalpies of liquid water and air water gas phases The mixing rule that used for the density of a liquid phase is 1 wI 5 4 2 22 Pa pa where p are partial densities The density for the gas phase is calculated from the ideal gas law with a correction for water vapor based on the steam
6. Relative tolerances for time step size control See tolerdt The tolerance for a particular variable is based on multiplying the relative tolerance value by the maximum magnitude of the variable over the entire problem domain Default depends on the module tolerconv lt var gt lt num gt Absolute tolerances for Newton Raphson convergence The Newton Raphson iterations is assumed to converge if the absolute magnitude of the change in primary variables between the previous and current iterations is less than the tolerance This tolerance is based on the values set by tolerconv and reltolerconv lxii APPENDIX These parameters are ignored if the rms NR conv option is set to on Default depends on the module reltolerconv lt var gt lt num gt Relative tolerances for Newton Raphson convergence See tolerconv The tolerance for a particular variable is based on multiplying the relative tolerance value by the maximum magnitude of the variable over the entire problem domain Default depends on the module cutbackmax lt num gt Maximum times that a time step is cutback because of lack of convergence of either Newton Raphson or preconditioned conjugate gradient iterations before the run is forced to terminate Default 30 rms time step lt on or off gt Turns on or off time step convergence based on root mean square of the estimate of the relative time truncation er ror If the error exceeds the tolerance set by rmstolde
7. It is often useful to be able to start a run from the results of another different run For instance one may want to use steady state conditions created from an initialization run as initial conditions for another run Another example is when modeling the effect of varying the duration of a remediation strat egy Restart files are written at various remediation times and subsequent no remediation runs are started at these different times to determine the effectiveness of remediation and the transport of remaining contaminants Restart records are written out to a file using the restart option in the output statement described below Note that records at different times can be written onto the same file The following statement is used to read in initial conditions from a restart file read restart file lt file name gt 4 USING THE COMPUTER PROGRAM lv time lt t0 gt read restart Here lt file name gt is the name of the restart file The program will pass restart records that corresponds to times that are less than lt tO gt It will read the first record whose time is greater than or equal to lt tO gt The run will start at simulation time equal to the value lt t0O gt unless the time statement is present in which case the initial time will start from the value set there It is an error to have both read restart and state statements Initial conditions read from the restart file can be overwritten for selected element
8. in the saturated zone as well as for the vadose zone If the water table of the saturated zone has a gradient a model with dimensions greater than one will obviously have to be created with the lateral boundary placed sufficiently far away where pressure heads which will be maintained consistent with the gradient using auxiliary elements If transport in the saturated zone does not need to be modeled such as in a one dimensional vadose zone model the bottom boundary at the water table can be kept fixed by an auxiliary element with aqueous phase saturation set to unity and pressure aqueous phase pressure set to the hydrostatic pressure equal to the half height of the element so that the water table will be at the interface between the auxiliary element and the vadose zone element immediately above it Concen tration can be set to some value often zero For a boundary condition of the third type the tortuosity factor can be adjusted so that the effective diffusion coefficient is equal to that of lateral dispersion in the saturated zone For a problem using the usnt module with a gas phase one has to be careful because an inconsistency may result from first or third type boundary conditions when there are multiple phases In particular it is possible for the gas phase to anomalously flow from the vadose zone into the saturated zone boundary element because the pressure in the element is not allowed to equilibriate with the gravitational head o
9. lt t2 gt List of output times file ext res File name extension of of restart output file restart forcetimes Force the model to hit the list of output times performs no output action Couttimes lt t1 gt lt t2 gt List of output times forcetimes output The contour option will write its output to a file with file name extension nvc The options elem history connect history and bc history are for writing the time history of an element variable a connection variable and boundary flux variables respectively and will all write to a single file with extension nvh Both the nvc and nvh files are read by the nview program a b c The list of output times must be strictly increasing order and the last output time must be less than or equal to the ending time of the simulation Using an astericks in the list of output times i e denotes that output will be at every time step The output file names are created by taking the part of the NUFT input file name up to the character before the first period and appending the characters specified by the file extension For example if the input file was called runi inp the NPREP input files would be run1 nvc and runi nvh and the restart file would be runi res Instead of specifying the file extension for the restart one can spec ify the full name For example one can use the statement file outtimes 4 USING THE COMPUTER PROGRAM
10. Is the problem isothermal or nonisothermal On the basis of the conceptual model a A NUFT module is then selected which can handle these processes efficiently The following criteria may be used to select the most efficient module a Does the problem involve nonisothermal conditions If yes use the usnt module b Is the problem saturated If yes use the ucsat module for flow and the usic module for contaminant transport c Is the problem unsaturated but with transport occurring only in the liquid phase If yes use the ustp module for flow and the usic module for transport d For all other problems use the usnt module for both flow and trans port Once we have identified the domain and its zones associated with dif ferent materials different initial conditions different distributed sources etc we establish the computational grid The latter consists of elements in rectangular or cylindrical coordinates according to the selected coordinate system The elements are formed by subdividing the x y and z axes in rectangular coordinates and r 0 and z in cylindrical coordinates The sub xxii APPENDIX divisions must be chosen such that element boundaries are on the boundaries of zones of homogeneous material properties boundary conditions homoge neous distributed sources and initial conditions As a rule consecutive subdivisions should not increase by a factor larger than two Boundaries of the zones may have to be rea
11. NUFT solves them In all cases it will be assumed that the solid matrix is station ary and nondeformable 0 0t 0 2 1 Saturated Flow The module ucsat is designed to solve both unconfined and confined prob lems of saturated flow The mass balance and flux equations for the flow of water are p k _Vy VP 2 1 AE Veq q z VP pgVz 2 1 iv APPENDIX Note that in NUFT p P For a compressible fluid the constitutive relationship for the fluid s den sity is p p P The fluid s coefficient of compressibility 8 is given by 1 Op 2 2 B SIP 2 2 Assuming o gt a Vp 2 3 r gt XV 2 equation 2 1 reduces to OP k OG Va asoy VP pgV2 2 4 Thus in the general case of flow of a compressible fluid in a nondeformable porous medium the variable is pressure P P x y z t For p constant 2 4 can be written in terms of the piezometric head h h y z t which serves as the variable h s V q So pgs q K VA 2 5 where S is the specific storativity Horizontal flow in confined and phreatic aquifers For essentially horizontal flow in a confined aquifer the variable is h h a y t The mass balance and flux equations are om V 0 Q 2 y t S S B T KB Q T VA 2 6 where B represents the aquifer s thickness and Qu represents symbolically the source term which usually takes the form of net recharge negative ext
12. be set to 0 6 In some cases value may have to be increased slightly if oscillations in concentrations develop Default 10 tolerdt lt var gt lt num gt Absolute tolerances for time step size con trol The primary variable name is set by lt var gt and the tolerance by lt num gt For example tolerdt P 1e3 C 0 01 C TCE 1e 6 S 0 05 S NAPL 0 01 T 0 5 tolerdt specifies that the pressure tolerance is 1e3 all mole fraction tolerances are 0 01 except for the component called TCE which is 1e 6 all sat uration tolerances are 0 05 except for the phase called NAPL which is 0 01 and temperature is 0 5 A time step is accepted if the absolute magnitude of the change in primary variables between the previous and current time step is less than the tolerance This tolerance is the larger of the absolute tolerances specified resulting from values given by tolerdt and reltolerdt The size of the time step of the next time step is based on estimating which value will cause the changes in primary variables a certain amount below the tolerance If the time step calculated at any Newton Raphson iteration is too large because the primary variables changed too much the time step size is cut back and the time step is redone If the rms time step option is set to on then the tolerance apply only to determining if a primary variable changes too much over the time step Default depends on the module reltolerdt lt var gt lt num gt
13. coordinates or cylind for cylindrical ones 2 The user must specify the location of the main coordinate system by choosing the coordinate origin and the direction of the coordinate axes with respect to the problem domain Typically the origin will be at a point of lowest elevation in the problem domain and the coordinate are chosen parallel to domain or subregion boundaries The components of the downward pointing vector lt down_x gt lt down_y gt lt down_z gt must be defined with respect to an auxiliary rectangular coordinate system For rectangular coordinates this auxiliary system is the same as the main system For cylindrical coordinates the x and y axes are defined as being perpendicular to the z axis and lying on the 0 0 and 0 90 deg planes respectively The down pointing vector need not be a unit vector because the model will internally normalize the vector If all of the components are zero then the model will set the gravitational terms within the balance equations to zero For the usual case the input item defining the direction of gravity is 0 0 0 0 1 0 3 The user must now define the computational grid The default shape of an element in NUFT in rectangular coordinates is a rectangle in 2 d do mains and a parallelepiped boxes 3 d domains In cylindrical coordinates it is an annular box This means that in all cases the edges of an element are along lines of constant coordinate values XXXi1 APP
14. flow factor lt tO gt lt f0 gt lt t1 gt lt f1 gt Optional specifies 3 table of time dependent flux factors 3 for this interval 3 default flux factor is one for all times can be used to shut off or open an interval or xlviii APPENDIX modify its flow characteristics setint End of specifications for this interval setint Specify next interval setint Can add more intervals intervals lt well name gt lt well name gt Specify next well lt well name gt wae 3 Can add more wells setwells Note the wells uses the permeability or hydraulic conductivity Kx of the completed element to compute its flow characteristics with the formation In particular for usnt NUFT uses the formula for well flow at a particular interval K kra da 2th Palpha well Pa elem PoG Zwell Zelem gt 4 4 Ha In te Tw where the subscript well and elem refer to the the center of the well interval and the centroid of the completed element respectively Here A is the interval length The permeability is set equal to Kx of the element A specified pressure production well lt well type gt is set to pres prod has the following parameters Table of time dependent pressure specified at the deepest interval of this well ptab lt t0 gt lt pO gt lt t1 gt lt p1 gt Optional friction loss factor pressure drop in well bore equals dp dl rho fac v v wher
15. interface between elements 1 and 2 The vectors t and n denote respectively the tangential and normal unit vectors to the boundary between the two ele ments This method has the disadvantage of excessive lateral dispersion at the contact between materials of highly dissimilar permeabilities with a large contrast in flow velocities The reason for this error comes from the as sumption that the dispersion coefficient varies smoothly between elements 1 and 2 whereas in reality it may have a discontinuous jump at their common boundary The method used in NUFT avoids this difficulty by using a treatment that allows the dispersion tensor to be discontinuous at element boundaries Let C be the concentration at the boundary AB between elements 1 and 2 The approximate dispersive flux through AB from element 1 to the boundary is disap Kni Cx C1 Ka Ca Cp 3 13 The flux from the boundary AB into element 2 is ipo Kn2 C2 Cx Ko Ca CB 3 14 where Kni n D n Lj i 1 2 3 15 Kiu n D t Lap i 1 2 3 16 We assume that to a good approximation these two fluxes are equal and therefore d d d AB TAB 2 A2 3 17 From these equalities we obtain an expression for C Substituting this expression into the expressions for Gap and qlz we obtain qiz Kr1Kn2 C2 C1 Kakn2 KniKi2 Ca CB Kni Kn2 3 18 This flux may be viewed as a generalization of a
16. is a list of some options that affect the numerical solution of the model and their default values time method lt method name gt Time step method Name of method can be fully explicit or fully implicit Default fully explicit for us1c and fully implicit for all other models time weighting lt num gt Linear weighting factor a defined in Subsection 3 1 Valid only if time step method is fully implicit Recom mended value for transport problems is a 0 6 which gives slightly less accuracy than a 0 5 Crank Nicholson method but is more nu merically stable Recommended value for flow problems especially for simulations used to obtain steady state conditions is a 1 0 Default 1 0 C itermax lt num gt Maximum number of Newton Raphson iterations before a time step is reduced and the iterations restarted Default 8 C iterbreak lt num gt Maximum number of Newton Raphson iterations be fore convergence is assumed Should be set to one if the equations that are being solved is linear because in that case only one iteration is necessary Default set to a very large number so that convergence is never assumed based on number of iterations but on other criteria linear solver lt option gt Method used for solving the system of linear equations during each step of the Newton Raphson method Some of the valid options are lublkbnd which is block banded LU decompo sition for a constant band width matrix d4vband which is the dire
17. lines must be drawn such that each homogeneous zone is defined as the set union of subsets of element groups c To specify boundary conditions the problem domain is extended by adding a strip or slab in three dimensions of boundary elements of constant width immediately outside the domain along every bound ary on which conditions of the 1st and 3rd kinds prevail No extension is required along no flow boundaries The width of the strip or slab can be taken to be the same as the size of the elements on the interior side of the boundary or any other size depending on the boundary condition to be implemented Zones of homogeneous material type are defined through the parameters in the mat mat statement The name of each zone is set by the lt zone gt parameter Note that names of elements within the zone will be of the form lt zone gt lt i gt lt j gt lt k gt where lt i gt lt j gt lt k gt are the indices i 1 nx j 1 ny k 1 nz 4 USING THE COMPUTER PROGRAM xxxiii of the subdivisions along the coordinate axes of the elements The name of the material within a zone is set by the lt mat gt parameter The following parameters lt imin gt lt imax gt lt jmin gt lt jmax gt lt kmin gt lt kmax gt indicate that the the zone contains elements with indices from 7 lt imin gt to lt imax gt from j lt jmin gt to lt jmax gt and from k lt kmin gt to lt kmax gt In rectangular coor
18. set in genmsh in order to set the flow lengths from the elements to the rest of model to zero In this way if one considers the elements as auxiliary elements one can place the flux at the boundary of the domain consisting of non auxiliary elements The following items apply a Linear interpolation is time is used in all tables b Positive flux is flux into the elements and negative flux is flux out of the elements c Fluxes have units of kg s except for usip and ucsat which have units of m s d Specific enthalpy is in units of J kg e Concentrations variable name is C lt comp gt are in mole fractions un less the statement input mass fraction on is present in which case the concentrations variable name is X lt comp gt are in mass frac tions The entire setcomp statement in phaseflux can be replaced by setcomp internal in which case concentrations and specific enthalpy will be those of each corresponding element set in range The setcomp internal option is the usual one if the fluxes are going out of an element i e are negative g Possible optional allocation methods are a Allocate the flux by volume allocate by volume xlvi APPENDIX The flux at element m equals the specified flux multiplied by the factor Um Xm Um where m runs over all elements specified in the range of the particular source term and Um are element volumes Allocate flux by flow area allocate by area l
19. stokes equations Comp Meth in App Mech and Eng 32 199 259 1982 Cuthill E and J McKee Reducing the bandwidth of sparse symmetric matrices ACM Proc of 24th National Conference New York 1969 Lee K H A Kulshrestha J J Nitao Interim Report on Verification and Benchmark Testing of the NUFT Computer Code Lawrence Livermore Laboratory Report UCRL ID 113521 1993 Nitao J J and J Bear Potentials and Their Role in Transport in Porous Media Water Resources Research to be published 1996 Nitao J J Reference Manual for the NUFT Flow and Transport Code Version 1 0 Lawrence Livermore National Laboratory Rep UCRL ID 113520 1996 Smolarkiewicz P K A fully multidimensional positive definite advection transport algorithm with small implicit diffusion J Comput Phys 54 325 362 1984
20. surface The tortuosity factor of the elements is set in rocktab to be equal to the ratio of the effective diffusion coefficient eddy diffusivity to the free diffusion coefficient The eddy diffusivity across the boundary layer is a function of wind velocity and vegetative cover the reader is referred to any mi crometerology text for details The capillary pressure of the elements must be set to identically zero The relative permeability of the gas 4 USING THE COMPUTER PROGRAM xli and liquid phases are set to one and zero respectively The porosity of the material is taken to be unity Ponded contaminant source of constant depth To model a ponded source of contaminated water with usip usic or usnt one uses a row of auxillary elements over the lateral extent of the desired pond area The permeability or hydraulic conductivity is set to be two orders of magnitude larger than that of the maximum permeability of the soil surface The porosity of the material is taken to be unity The relative permeability of the liquid phase is set to one and a linear relative permeability k g S is used for the gas phase Any normalized Kg values are set to zero The aqueous phase satu ration of the auxiliary elements is kept equal to one For the usnt module the absolute pressure P for usnt is set to the equivalent pres sure of a head of water plus the atmospheric pressure say 1e5 Pa For the ustp module the pressure head H is set to the depth
21. t1 gt lt v1 gt phasefactor 4 USING THE COMPUTER PROGRAM XXXIX lt bc name gt End of this boundary condition lt bc name gt Name of boundary range lt range gt lt range gt Range of elements clamped Specifies that primary variables 3 are set to those set by initial condition using tt state or tt read restart options factor factor Optional component flux factors phasefactor phasefactor Optional phase flux factors lt bc name gt One can add more boundary 33 conditions here bctab End specification of boundary conditions The boundary condition name lt bc name gt is unique and is used for output purposes such as for keeping track of the total flux across the boundary segment The range statement specifies the range of auxiliary or bound ary elements for the boundary segment The basephase statement specifies the phase for which the concentrations will specified in the tables state ment The tables statement gives the values of the primary variables in the auxiliary elements in the range statement Values are given in the form of tables at times lt tO gt lt t1 gt and values lt valueO gt lt valuel gt at these times Linear interpolation is used to calculate values between the time entries in the table The times must be in strictly increasing order with the last time greater than or equal to the ending time as given by th
22. tables The enthalpy of a phase is approximated by the mixing rule ha X wih 2 23 0 3 Numerical Model 3 1 Finite difference discretization of partial differential equa tions In both the ucsat and usip modules NUFT numerically solves a single partial differential equation the mass balance equation for water The usic module solves a single mass balance equation for the transport of a single component In the usnt module NUFT solves a coupled set of balance equations one for the mass of each transported component For the non isothermal case an energy balance equation is also included in the coupled set As we have seen above each balance equation has more or less the general form abe ot where 6 S is the volumetric fraction of the fluid e is the density of the considered extensive quantity J is the non advective flux of the considered V be V J 0T 3 1 extensive quantity and I describes the strength of the source as quantity per unit volume of phase For a fluid s mass e is replaced by the mass density p In order to discretize this partial differential equation in time and space the time domain is subdivided into a sequence of time steps t n 0 1 and the space domain is partitioned into finite difference cells or elements Specifying an initial time step the length of subsequent x APPENDIX time steps is automatically controlled by the program using user specified tolerances f
23. used by the ORTHOMIN algorithm Recommended value is 25 for d4 and ilu and 30 for bgs The tolerance for convergence is set by the toler param eter Recommended value is 0 001 for most problems In some cases a smaller value 0 00001 is needed if frequent instances of Newton Raphson nonconvergence is observed The maximum number of al lowed iterations is set by the itermax parameter If this number of iterations is reached the model cuts back the time step and retries the entire Newton Raphson method Recommended value is 200 for d4 method and 500 for ilu and bgs methods Default none because pcg is not default for linear solver ilu degree lt num gt The degree of fill in in the incomplete LU decompo sition algorithm relevant to the ilu and d4 preconditioning options Recommended value is 0 for saturated flow and transport problems 1 for nonlinear problems including when running the usip module and for most problems using the usnt module and 2 or more for severely nonlinear problems including many multiphase problems us 4 USING THE COMPUTER PROGRAM lxi ing the usnt module Default 1 artificial diffusion lt off or on gt Turns on or off the artificial dif fusion method of Brooks and Hughes 1982 for reducing numerical dispersion Valid options are on or off Default off artificial diffusion coef lt num gt Sets artificial diffusion coefficient Used for transport calcuations Recommended that time weighting parameter
24. v refers to v th iteration The symbol J is the Jacobian matrix of R The initial iterate v is taken from the solution at the previous time step Note that we must solve a system of linear equations at each iteration NUFT has two different options for determining when convergence of the iterations has been sufficiently reached One method is to stop when changes in solution between two iterations are sufficiently small The other method is to stop when the current root mean square of the residual R is smaller than a specified tolerance 3 2 Switching of primary variables The set of primary variables that is solved in a particular element depends on the particular phases that are present or absent in that element For ucsat the primary variable for a saturated element is the pressure head of water with atmospheric pressure used as zero datum For an un saturated element the primary variable is the fractional portion of the void space occupied by water within the element We refer to this fraction as element saturation e An unsaturated element is switched to saturated in the course of a simu lation if the element s saturation becomes greater than or equal to unity e A saturated element is switched to unsaturated if the pressure head in i e at the centroid of the element becomes less than the element s vertical thickness For ustp the primary variables are pressure head with zero datum at at mospheric pressure and satura
25. word with as its first character or the end of file can be used as the terminator instead of end Multiple variables can be present in a single file the model will search through file for the first occurrence of the primary variable name before reading the element data that immediately follows Examples state C contam by ijk file initc dat C water 1 0 atm 0 02 C contam by ztable range WT 3 height 33 above 3 model base C contam 0 0 0 001 18 0 0 01 20 0 0 01 21 0 0 0 S liquid by key 1 0 liv APPENDIX S liquid by ztable range VZ 3 height 33 above 35 WT 12 12 12 13 20 liquid S 0 9 5 4 3 Oo OWN OO oo 0 0 F P by key 1 e5 In this fictitious example the concentrations C contam of a contaminant are first read from a file initc dat In the next specification the concentra tion of the water component C water is set to 0 02 for elements that start with the characters atm and set to 1 for all other elements denoting the fact that the concentration is calculated internally Then the concentrations for the contaminant at elements that begin with WT are overwritten by values calculated from a table Then the liquid saturation liquid is set to 1 0 for all elements The liquid saturation for those elements that start with VZ are then overwritten by a table Finally the pressure is to 1le5 for all elements 4 13 Using the restart option
26. 33 z direction mat Begin assigning zones and materials 33 Set lt zone gt to have material lt mat gt lt zone gt lt mat gt lt imin gt lt imax gt lt jmin gt lt jmax gt lt kmin gt lt kmax gt Repeat for as many different materials are present lt zone gt lt mat gt lt imin gt lt imax gt lt jmin gt lt jmax gt lt kmin gt lt kmax gt mat End zones and materials null blocks Specify null or inactive elements lt imin gt lt imax gt lt jmin gt lt jmax gt lt kmin gt lt kmax gt lt imin gt lt imax gt lt jmin gt lt jmax gt lt kmin gt lt kmax gt 4 USING THE COMPUTER PROGRAM XXX null blocks 3 End null blocks 3 SPECIFY BOUNDARY ELEMENTS AND CONDITIONS In what follows replace lt range gt by the name assigned to a boundary segment bc 1st type range lt range gt lt range gt bc 1st type bc 2nd type range lt range gt lt range gt bc 2nd type bc 3rd type range lt range gt lt range gt bc 3rd type Following is optional If present then the hydraulic conductivities or permeabilities in the x y z directions 33 given in rocktab are used If not present the value 33 in the x direction is used for all three directions anisotropic genmsh End mesh specification Let us explain 1 The user must first select the type of coordinate system by setting coord type to either rect for rectangular
27. APPENDIX by John J Nitao Lawrence Livermore National laboratory Livermore CA 94551 USA Many flow and transport problems of practical interest cannot be solved analytically They must be solved numerically on a digital computer using a computer code or program Ideally for the purpose of solving the forecasting models discussed in this book a computer code should be able to solve a wide range of models in both the unsaturated zone and the saturated one including nonisothermal as well as isothermal cases The obvious reason for this statement is that as emphasized throughout the book we have essentially only one partial differential equation that appears in all models viz the equation that describes the balance of an extensive quantity Also the boundary conditions are essentially the same for all extensive quantities In most cases a model would involve the transport of more than a single extensive quantity so that a code should be able to solve a number of p d e s simultaneously However the actual numerical simultaneous solution of a number of p d e s and the transformation of such solution into a practical and useful computer program may require many thousands of code lines NUFT Nonisothermal Unsaturated saturated Flow and Transport serves as an example of a code that can be applied to a large number of flow and contaminant transport problems We have included this code with the book in order to facilitate the soluti
28. C NZP concentrations The disappearance of a phase is triggered by a saturation as a variable becoming less than or equal to zero Appearance of a phase is flagged when the sum of trial mole fractions of a phase is greater than or equal to unity during a Newton Raphson iteration Note that in order to generate the derivatives for the Jacobian matrix in the Newton Raphson method the model must consider the appropriate set of primary variables at each element Example of primary variables for an air water system Consider a system consisting of air a and water w as the only components Of course air is a mixture of many components but here we assume that the problem being modeled allows treating the individual constituents as indistinguishable so that air can be considered as a single component The gaseous phase g is a mixture of air and water vapor The aqueous phase is a mixture of liquid water and dissolved air The number of primary variables is equal to NC 1 3 Examples of primary variables are e If S 0 and Se 1 NZP 1 then we use py wy and T e If S 1 and S 0 NZP 1 then we use py w and T e If S Se gt 0 NZP 2 then we use py Se and T Note that in the last case the w s are defined by the thermodynamic equi librium assumption and are functions of p and T through the equilibrium partitioning coefficients xiv APPENDIX Figure 3 1 Definition sketch for numerical dispersion Exam
29. ENDIX The grid is defined by specifying the size of the elements in each coor dinate direction through the parameters lt dx_1 gt lt dx_2 gt lt dx_nx gt lt dy_1 gt lt dy_2 gt lt dy_ny gt lt dz_1 gt lt dz_2 gt lt dzmnz gt where nx ny nz are the number of elements in each coordinate direction The subdivisions start at the origin of the coordinate system For example dx 1 100 10 50 5 30 1 50 dy 10 20 1 40 dz 1 100 means that along the x axis we have 17 elements of the indicated sizes one of Ax 100 m then 10 times Av 50 m finally one of Az 50 m Along the y axis we have 11 elements 10 times Ay 20 m then one Ay 40 m Along the z axis we have only one element Az 100 m The z coordinate of the element centroids are x lt dx_1 gt 2 rg ay lt dx_2 gt 2 For cylindrical coordinates the coordinates x and y corresponds to the r and coordinates respectively with the angle in units of degrees The following steps describes how to define a grid for a two dimensional problem in rectangular coordinates However this can be extended to three dimensional problems or to cylindrical coordinates in the obvious way a Draw the problem domain and enclose it in the smallest possible rect angle with sides that are parallel to the selected x and y axes b Subdivide the rectangle into elements by lines parallel to the z and y directions that extend throughout the rectangle The
30. GRAM xxxvii 4 8 Specifying component and phase properties The usic and usnt modules requires that component property parameters e set For ustc only the free binary diffusion coefficient m s is set by the statement compprop lt comp gt lt phase gt freeDiffusivity lt value gt compprops where lt comp gt and lt phase gt are the names of the component and phase respectively set in the init eqts statement For usnt the following parameters are set compprop lt component gt Name of component set of parameters 3 is needed for each mass component intrinsic Parameters that do not depend on which phase go here MoleWt lt value gt Molecular mass g mole intrinsic lt phase gt Properties of component that depend on which phase it is in go here Keq Equilibrium partitioning coefficient freeDiffusivity Free diffusion coefficient rhoc Partial mass density Centhc Partial specific enthalpy not present for isothermal problem lt phase gt lt component gt compprop We refer the reader to the NUFT user s manual for the usnt module for details For usnt the following phase dependent parameters are set phaseprop lt phase gt Name of phase set of parameters needed for each fluid phase rhoP Mass density of the phase viscosity Fluid viscosity of the phase CenthP Specific enthalpy of the phase n
31. HE COMPUTER PROGRAM XX X lt flow module name gt Name of flow module e g ucsat usip or usnt Flow module input parameters go here lt flow module name gt The following item is included in a transport module lt transport module name gt Name of transport model e g usic Transport module input parameters go here lt transport module name gt 4 5 Specifying the names of phases and components In the ucsat and usip modules liquid is always the name of the phase and water is the name of the component For the usic and usnt modules the phase and component names must be specified by the following statements init eqts Define components and and component names components lt compO gt lt comp1 gt Names of components phases lt phase0 gt lt phase1 gt Names of phases wetting phase lt phase gt The wetting phase partitioning allowed between this phase and 33 the solid primary phase lt phase gt The primary phase other 3 phases are computed w r t this phase using capillary pressure thermal Present only if problem 33 is not isothermal isothermal Present only if problem 33 is isothermal init eqts Only the components and phases parameters are set for ustc module The following points apply to the usnt module a Either the isothermal or thermal must be set depending on whether the problem is assumed to be isothermal or not An
32. Using the bc 3rd type command NUFT optionally sets B 1 so that any element width may be used In this case c 1 K i e the resistance will be the reciprocal of the hydraulic conductivity For multiphase problems we have that the flux at the boundary of the a phase is dom A E pa Pao Pagle 2 3 22 Qa where we specify the transfer coefficient A by the value of permeability K The value of k is based on upstream weighted values as defined by 3 4 Note that if the saturation of a phase in the auxiliary element is zero then there will be no flow coming into the boundary of that phase because kpa is zero A boundary condition of the third type occurs in a contaminant transport problem when we assume the presence of a well mixed zone on the external side of a boundary This kind of boundary is discussed in Subs NUFT uses condition rewritten here for convenience in the form c qarn A D Ven 3 23 in which c denotes the concentration in the external domain n denotes the normal to the boundary and A plays the role of a transfer coefficient In the numerical solution NUFT approximates 3 23 by the equation c c q n A Darc e 3 24 where the value of Dp is the harmonic average between the values of the auxiliary and adjacent interior elements are used The value of D can be adjusted by changing the the tortuosity factor of the material type of
33. XIV APPENDIX The null blocks command specifies which elements are inactive The bc 1st type bc 2nd type and bc 3rd type commands are used for the auxiliary elements the added strip of elements just outside the prob lem domain They indicate which type of boundary condition should be implemented along a specified boundary segment In numerical implementation of the boundary conditions the flow length from the center of elements assigned in bc 1st type or bc 2nd type to the boundary i e any neighboring element that is not an auxiliary element is set to zero The flow length in elements indicated in the bc 3rd type command is unity If elements in the bc 1st type and bc 3rd type commands are not specified also in the bctab command to be described later or if elements in the bce 2nd type command are not specified in the srctab command the program will indicate a fatal program error 4 7 Material property parameters The parameters that describe the properties of each material type are given in the statements rocktab lt material name gt parameter values lt material name gt lt material name gt parameter values lt material name gt rocktab Here lt material name gt is the unique name given by the user to a given material type The input parameter values for each material type depend on the particular module We will briefly describe the parameters here The user is referred to the appropriat
34. action of fluids These elements are similar to those used for a ponded source except that in most cases the inflow elements should have the tortuosity factor set to zero so that there will be no binary diffusion flux coming out of them Otherwise suppose for ex ample that the concentration of a contaminant is specified to be zero and that there is an inflow advective flux of contaminant coming into the element The very sharp gradient in contaminant concentration that results between the soil element and the inflow element will often cause numerical instabilities and oscillations Downstream boundary condition A common problem that is modeled is one dimensional transport of a contaminant say through 4 USING THE COMPUTER PROGRAM xliii a soil column An element at the inflow side is kept at constant pres sure and concentration The width of the element can be made very small or the bc 1st type option can be used At the outflow side the element is kept at constant pressure and zero concentration but with zero tortuosity factor for the same reason as stated for the well inflow boundary condition Water table In building a vadose zone model one may or may not wish to include the saturated zone depending on the domain of interest If the saturated zone is included rows of elements at the lower part of the domain is placed with initial aqueous saturation equal to one An initialization run may have to be made to equilibriate the pressures
35. also possible A flux can also be applied inside the considered domain by specifying a flux into an element at its centroid NUFT allows the user to specify either the mass flow rate of a particular component or phase In order to specify a flux of a phase the concentration of the phase must also be specified in addition to the rate of flow of the phase When withdrawing fluid from an element the user can request that the concentration of the phase is that of the producing element The boundary condition of the third type in a flow problem occurs when ever we have a semipervious thin domain between an external domain of specified head or pressure and the considered domain In other words a xviii APPENDIX resistance exists to the flow from the external domain in which the head or pressure is specified into the considered domain For saturated flow refer ring to 3 21 or rewritten here for convenience in terms of h as the variable qn Z A h ho 3 21 Cr where is a transfer coefficient and c is the hydraulic resistance of the semipervious membrane In NUFT the resistance c is made equal to B Kk where B is half the width of the auxiliary element distance from centroid of auxiliary element to the boundary and K is the hydraulic conductivity of the material assigned to that element By selecting the hydraulic conduc tivity for a given element width the desired value of c can be achieved
36. at are not just a single number but a list of numbers or a list of other associated input parameters For example sand K 1 2e 12 porosity 0 2 sand assigns the properties hydraulic conductivity K and porosity for a material type called sand Parenthesis must always come as a pair they play an important role in the input file NUFT will check and issue a diagnostic message if either or are missing The second form shown above is therefore recommended for complicated statements that contained nested parentheses because NUFT will then indicate the exact location of any mismatched parentheses in the issued diagnostic message We shall use the second form in the examples below wherever a statement contains a nested set of parentheses An input parameter may specify types of correlations For example kr krVanGen m 0 2 alpha 1e 4 Sr 0 2 kr specifies the set of parameters to be used for the relative permeability kr 4 USING THE COMPUTER PROGRAM XXV as defined by the Van Genuchten correlation Following are some rules and conventions used in a NUFT input file a Spaces tabs and line breaks are ignored within an input file except to separate data or key words Exponents in floating point numeric values can be in e e g 3 4e3 or E e g 3 4E3 format but d or D FORTRAN style format is forbidden The ordering of input state ments within another statement is irrelevant unless
37. auxiliary elements next to the external side of the considered boundary As we shall see the type and material of the auxiliary elements depends on the boundary condition type The boundary condition of the first type is implemented in NUFT by specifying the possibly time varying primary variables of the auxiliary ele ments next to the considered boundary segment Although these primary variables are not calculated during the simulation they take part in the calculation of the fluxes through the boundary The values of the primary variables can either be specified as a piecewise linear continuous function of time or by the clamped command the primary variables will be equal to those specified in the input file as initial boundary values The type of porous medium material in an auxiliary boundary element is usually set 3 NUMERICAL MODEL xvii to be the same as in the adjacent internal element The range of elements in the bc 1st type command described below must include the auxiliary elements in order to internally make the flow length from the center of these elements to the boundary equal to zero This action forces the primary vari ables to be set at the boundary edge by effectively moving the centroid of the auxiliary element to the midpoint of the common boundary surface between the auxiliary element and the adjacent interior element The width of the auxiliary elements as well as their volumes are immaterial NUFT has the option o
38. ch is specified their minimum and maximum 2 J k indices i e i lt imin gt lt imax gt j lt jmin gt lt jmax gt k lt kmin gt lt kmax gt In the documentation of the NUFT input that follows the index range method can be used anywhere that it is stated that the range method is used h The following statement include lt file name gt will paste the contents of a file in the current directory into the input file It can be placed anywhere in the input file where a complete statement with matching parentheses is expected i The statement define lt symbol gt will cause the stated symbol to be defined It is used in conjunction with the ifdef and ifndef statements in order to conditionally activate selected sections of the input file That is if the stated symbol with name lt symbol gt is defined then the input statements within ifdef lt symbol gt 4 USING THE COMPUTER PROGRAM xxvii ifdef will be read as input data If lt symbol gt is not defined by a define statement then the statements will be ignored Conversely the state ments inside of the following ifndef lt symbol gt ifndef will be read if and only if the symbol with name lt symbol gt is not defined A special symbols is the one called INIT It is considered defined if and only if the program nuftini is used to run the input file The program called nuft does not define the symbol INIT This fact is the onl
39. cities are given by darcy s law ka Va Tin VPa pogV2 The phase pressures are given by the capillary pressure relationships For example for a system with gas g and liquid phases Pe Pg Pel Se T 2 15 For a NAPL n and aqueous phase system Pe Pn Pe Se T 2 16 assuming that the porous medium is water wet For a gas NAPL aqueous phase system Pe Pn Peh Sag T Pn Pg Pel Sg T 2 17 To determine all w s we require partitioning coefficients that relate component concentrations between each pair of phases ng Kien 2 18 vill APPENDIX In using this relationship note that the mole fractions n2 nh must be con verted to mass fractions wl w3 NUFT has hard wired partitioning coefficients for the partitioning of water between the gas and aqueous phases based on at Dsat T pg exp Y M RTabspe 2 19 where w is the matric suction potential for the aqueous phase Tabs is absolute temperature and Psat is the saturation pressure from the steam tables The exponential term is often called vapor pressure lowering and is derived from the porous media generalization of Kelvin s law Nitao and Bear 1996 For TCE NUFT also has hard wired values of Henry s coefficients for the partitioning between the gas and aqueous phases A general functional form for partitioning coefficients is KY A py B C pg exp D Tavs
40. ct solver with D4 ordering Behie and Forsyth 1984 pcg which is the preconditioned conjugate gradient method with various precondition Ix APPENDIX ing options The option lublkbnd is best for small problems d4vband for intermediate problems and pcg for large problems especially in three dimensions Default lublkbnd pcg parameters Parameters for preconditioned conjugate gradient method Following is the complete format pcg parameters precond lt option gt north lt num gt toler lt num gt itermax lt num gt pcg parameters Some valid values for precond are bgs which is preconditioning by a single iteration of the symmetric block gauss seidel method ilu which is the incomplete LU decomposition method d4 which is the incomplete LU decomposition method applied to a system that has been reduced by D4 ordering and comb which is the combinative method The reader is referred to Behie and Forsyth 1984 for a description of these methods and the ORTHOMIN algorithm which is version of the conjugate gradient method for nonsymmetric matri ces that used by the model The bgs method is recommended for relatively easy problems without severe nonlinearities such has satu rated flow and transport problems The d4 method is recommended for most other problems The ilu method uses less memory then d4 but usually requires more conjugate gradient iterations to converge The north parameter is the number of orthogonalizations
41. cular there are two types of graphics output files XXIV APPENDIX The file ending with nvc standing for NUFT variables contours contains information for generating the data required for mapping contours say of H S P while the file ending with nvh standing for NUFT variables his togram contains the information required for plotting requested histograms of variables say H t S t at specified locations Step 3 To obtain graphical output whether contours or history plots we use the program nview which is executed by clicking on its icon In the file menu we select the name of the input file for the run without the file s extension For example selecting nview p911 will search the current directory for the files p911 nve and p911 nvh and read them in accordance with the user s specifications concerning the files to be plotted 4 3 Input syntax and format We shall now describe the key parameters in the NUFT input file for those who wish to use a text editor in order to create or edit an input file The general syntax for an input item in the input file is either lt key word gt lt values gt or lt Key word gt lt values gt lt key word gt where lt key word gt is a key word for an input parameter and lt values gt are the values assigned to that parameter For example porosity 0 2 would set an input parameter that is called porosity to a value of 0 2 An input parameter may have values th
42. ded in Chap are provided with the NUFT software distributed with the book Each such file can be modified by using any text editor to create a new input file for any other set of parameters as long as the problem belongs to the same class of problems viz has the same structure of input file For example if we have a problem of saturated flow in a confined aquifer we may modify the file p911 inp provided with the software to be used for solving the problem described in Subs and name the newly created file newname inp e Use the NUFT preprocessor nprep by simply clicking on its icon which 4 USING THE COMPUTER PROGRAM xxiii should have been created after its installation The program will then display a menu from which the user first selects the desired problem say P9xy He may then modify parameters on the screen and save the file un der the same name P9xy inp or under any desired name newname inp By using the command nprep the user may later reopen the saved file and make further modifications or create a new file While preparing the input file the user is referred to the documentation within the program for help or further information on how to prepare nprep NUFT is highly flexible and has many more input options than described here More details may be found in the NUFT Reference Manual Nitao 1996 which is included with the distributed code Step 2 Once an input file filename inp has been p
43. dinates the elements so defined comprise a rectangular zone however more complicated zones may also be defined Multiple zone definition statements with the same zone name and or material type are permitted with the elements within any zone overriding effect of any previous statement of zone definition Thus two options are possible for the definition of zones Option 1 A zone is defined as a union of rectangular zones Option 2 A zone is defined by overlaying it i e writing the statement defining it later in the file on top of another zone Note the equivalence of the following two examples According to Op tion 1 mat Ri sand 1 nx 1 3 1 nz Ri sandi 1 4 ny 1 nz Ri sand2 3 7 ny 1 nz Ri sand4 4 4 ny 1 nz Ri sand 5 nx 4 6 1 nz Ri sand 8 nx 7 ny 1 nz R2 silt 2 3 4 6 1 nz R3 clay 57 7 ny 1 nz mat According to Option 2 mat Ri sand 1 nx 1 ny 1 nz First set all elements as R1 R2 silt 2 3 4 6 1 nz Cut out R2 R3 clay 57 8 9 1 nz Define R3 as the union of R3 clay 67 7ny 1 nz two rectangular regions mat Obviously there is no uniqueness and other subdivisions are possible It is obvious that in the above example Option 2 is preferable Note that the model will accept the symbols nx ny and nz instead of the maximum element indices in the z y and z coordinate directions respectively Also simple arithmetic expressions such as ng 1 and ny nz 1 are also allowed XX
44. djusted to accommodate this requirement Boundary elements are an exception to the above rule unless the boundary condition is of the second kind The default boundary condition in NUFT is a no flow boundary Bound ary conditions of the first and third kinds are implemented by adding a row or column of auxiliary elements outside the domain next to boundaries In Subs 3 5 we have described how these auxiliary elements are used to im plement these two kinds of boundary conditions 4 2 Translating the conceptual model to an input file The conceptual model must be converted to input parameters that can be read by the NUFT program These input parameters are specified in the NUFT input file The following steps describe how an input file may be prepared Step 1 Prepare the input file which contains input parameters for the NUFT program These parameters include the geometry of the computa tional finite difference mesh material properties boundary conditions and initial conditions In each case the particular set of input parameters to be used depends on the model that is selected Nevertheless the key parameters are similar for all models Each problem has its own input file Given a problem say P9xy i e described in Subs 9 x y there are three options for creating the required input file e Use the information in the NUFT Reference Manual to create the required input file filename inp e Input files for the problems inclu
45. e NUFT user s manual for more details For the ucsat module the following parameters must be set for each material type porosity lt value gt Fractional porosity Kx lt value gt Ky lt value Kz lt value gt Sat hydraulic conductivity m s 3 in the x y and directions 33 The Kx value is used for all directions unless the anisotropic option is set in genmsh For the ustp module the parameters are porosity lt value gt Same as used in ucsat above Kx lt value gt Ky lt value Kz lt value gt Same as used in ucsat above 4 USING THE COMPUTER PROGRAM pe lt function gt lt parameters gt 33 kr lt function gt lt parameters gt 33 The parameters for the ustc module are porosity lt value gt 33 Kd lt value gt 33 tort 33 lt function gt lt parameters gt 23 tort 33 dispL lt num gt 33 gt dispT lt num gt XXXV Capillary pressure head function m Relative permeability function Same as used in ucsat above Normalized linear isotherm adsorption coefficent unitless equal to bulk soil density times the usual Kd divided by porosity equal to one minus the retardation factor Tortuosity factor unitless equal to reciprocal of tortuosity multiplies binary diffusion coefficient Longitudinal dispersivity m the dispersion on statement must also be present for dispersion Transverse dispersivity m F
46. e tstop record If a sharp change in primary variable is specified the forcetimes statement in the output record should be used so that the model will not take too large of a time step and jump over the sharp change Instead of specifying the primary variables by tables one can use the clamped statement which will tell the model to use the primary variables set by one of the initial condition commands This option is especially useful when reading in steady state initial conditions read from a restart file and one wants the boundary to be maintained at steady state conditions The optional factor statement tells the model to multiply the compo nent fluxes calculated between the specified boundary elements and the rest xl APPENDIX of the elements not including other boundary elements by time dependent factors calculated from the associated tables using linear interpolation in time Note that there has to be a table present for every component flux corresponding to a separate factor for each component The phase factor is similar to the factor statement except that it multiplies the phase fluxes between the boundary elements and the other elements There is a time dependent factor and its associated table for each phase If both phase factor and factor statements are present the net result will that a component flux in a particular phase will be multiplied by the factor from phase factor for that phase and by the factor from factor for the com
47. e v is flow velocity 3 default value is 0 0 fricfac lt fac gt 3 Optional table for computing time dependent factor multiplying all phases default is one at all times using this option an entire well can be shutin 4 USING THE COMPUTER PROGRAM xlix by setting the factor to be zero at a particular point in time well factor lt t0 gt lt f0 gt lt t1 gt lt f1 gt The following are input parameters for a specified pressure production well lt well type gt is set to pres inj pres inj 3 name of injected phase phase lt phase name gt Table of time dependent pressure of the 33 shallowest interval for this well ptab lt t0 gt lt pO gt lt t1 gt lt p1 gt Table for fluid mass density in well bore kg m3 density lt tO gt lt rho0 gt lt t1 gt lt rho1 gt Following table of specific enthalpy is for thermal models only J kg enthalpy lt tO gt lt HO gt lt t1 gt lt H1 gt Optional friction factor for well bore flow see above fricfac lt num gt Optional well factor see above well factor lt t0 gt lt f0 gt lt ti gt lt f1 gt Tables for time dependent concentrations of mass components 3 in injected stream units are in mole fractions unless input mass fractions on is present in which concentrations are in mole fraction lt comp name gt is name of 3 the particular component conc lt comp name gt lt t0
48. els that are described in the various sections of this book focusing on the prob lems presented in Chap No prior experience with numerical methods or the use of computer codes is required 1 General Description 1 1 Objectives NUFT was written in order to solve subsurface environmental problems especially problems that require the prediction of contaminant migration e g the analysis of remediation methods and the study of nuclear waste isolation problems The original emphasis was on the transport of multiple components in a multiphase flow system It has recently been extended to make it applicable to a wide variety of saturated and unsaturated flow and contaminant transport problems However the more experienced user will recognize how the general case is reduced to the particular cases of saturated flow single phase flow under saturated or unsaturated flow conditions etc Essentially NUFT solves the partial differential equations that describe the mass balances of multiple components transported within multiple phases that together occupy the the void space of a porous medium domain The components and the phases are assumed to be in local thermodynamic equi librium as described in by Bear and Nitao 1995 The results of the calcu 2 MATHEMATICAL MODEL iii lations using NUFT are either time histories of concentrations saturations and fluid pressures at different locations within the problem domain or spa tial di
49. es midpoint weighting for the concentration in the advective term together with an artificial diffusion term that is a function of the grid Peclet number Alexander and Hughes 1982 e Using the flux correction method suggested by Smolarkiewicz 1984 which uses upstream weighting together with an anti diffusive term that can cels a portion of the first order numerical dispersion flux 4 Using the Computer Program 4 1 Developing a Conceptual model of a problem Before preparing the input parameters for the code the conceptual model Subs of the considered problem must be developed We start by identifying the problem domain The domain must contain not only the region of direct interest but also any possible larger domain of influence for the time span of the simulation We identify subregions in which we may need a more detailed solution and others in which a less accurate solution will suffice Sometimes we have to take into account avail ability of information on boundary conditions when selecting the domain s boundaries Within the problem domain we identify regions within each of which the material porous medium is homogeneous with respect to relevant proper ties We then identify boundaries and and conditions on them We divide the boundaries into boundary segments on each of which conditions are homogeneous i e the same head the same flux etc prevails along the entire segment We identify internal
50. f gas in the vadose zone but is kept at a constant value Even if the pressure at the water is set correctly any changes in gas pressure such as from vapor extraction or barometric pressure changes would cause unphysical gas flow One way to remedy this problem is to use the factor option to turn advective gas flow off at the water table Another solution is to add at least one row of saturated elements above xliv APPENDIX the auxiliary row of elements In this case gas cannot flow into the water table elements because any changes in gas pressure results in slight changes in the vertical position of the water table equal to the change in gas pressure 4 10 Specifying boundary condition of second type and source terms The following statements specify sources both as 2nd kind boundary condi tions where a flux is specified along a boundary as point sources and sinks actually these sources and sinks are indicated for elements rather than points and distributed sources over a range of elements e g infiltration in a 2 d phreatic aquifer model srctab phaseflux Present if one wishes to specify a 33 a given flux of phase name lt source name gt 3 Name of source term used in flux output option phase lt phase name gt Name of phase for which flux is specified range lt range gt lt range gt Range of elements table lt t0O gt lt q0 gt lt t1 gt lt qi gt Table of flux of the phase
51. f turning on or off the flux of specified components and or phases across a first type boundary This option can be used for turning on or off a contaminant source at a boundary A boundary condition of the second type is also implemented through the use of auxiliary boundary elements or a strip of such elements on the external side of the boundary However in this case the primary variables of the auxiliary elements do participate in the calculations The specified flux through the boundary is specified as a rate of discharge into each auxiliary element Again the width of the auxiliary elements is immaterial since the user by using the command bc 2nd type causes NUFT internally to move the centroid of the auxiliary element to the boundary thereby causing the specified flux to along the common boundary surface between the auxiliary and adjacent interior element The pore space volume in these elements should not be much larger than that in the neighboring element Otherwise there may be a significant storage effect which will result in a time lag in the flux across the boundary The porosity can be adjusted if necessary to change the volume of the pore space For the modules us1p and ucsat the discharge rate through the bound ary is specified in terms of m element s For usnt the flux is specified in kg element s Both rates may be constant or time dependent An option to specify fluxes per unit flow area of boundary surface are
52. gt lt CO gt lt t1 gt lt C1 gt lt comp name gt lt t0 gt lt CO gt lt t1 gt lt C1 gt conc 4 12 Specifying initial conditions The initial conditions to the model can be specified by either of two methods a By specifying values using the state option b By reading in a restart file created from another simulation run This option is described later Only one of the above methods can be used We first describe the state option It is of the following form l APPENDIX state lt primary variable gt lt method gt lt specification gt lt primary variable gt lt primary variable gt lt method gt lt specification gt lt primary variable gt state Here lt primary variable gt is the name of each primary variable The pri mary variables of each module has already been given The name of the primary variable for the ucsat and usip modules are a H denotes the pressure head in meters of water relative to atmospheric conditions That is setting H to zero is the same as atmospheric pres sure head Optionally if the statement input by piezohead on is present at the main input level then the variable H is piezometric head b S1 denotes liquid saturation Since only one of the above variables can be the primary variable of a partic ular element some rules for deciding which variable is being set is required In the usip module the following rules apply in the following order of prior it
53. harmonically averaged flux 3 4 Treatment of a phreatic surface as a boundary NUFT employs the following approximation to advance the position of the phreatic surface from one time step to the next If h is the average phreatic xvi APPENDIX surface height in the mth element the discretization of 2 9 is Asm pess B hp APH ER Aom NA 3 19 where A is the area of phreatic surface in the element qem is the z component of the fluid s flux just below the phreatic surface Assuming incompressibility we have e So art 3 20 meENm where qm m is defined to be the fluxes out of the element into neighboring elements m Nm 3 5 Specifying boundary conditions Three types of boundary conditions may be applied at any given boundary or segment of a boundary e Boundary condition of the first type Dirichlet condition i e a boundary along which the value of a variable is specified at all times e Boundary condition of the second type Neumann or flux condition i e a boundary along which the flux normal to the boundary is specified at all times e Boundary condition of the third type Cauchy or mixed condition in which a combination of the variable and its gradient are specified How does NUFT implement these boundary conditions In NUFT the default boundary is one of no flux or impervious All other types of bound ary conditions are implemented through the addition by the user of a strip of
54. heat solid density lt value Grain density of solid kg m3 tcond Bulk thermal conductivity lt function gt lt parameters gt 33 CW m C tcond The following important notes are for the usnt module a In all of the above a property can be set to a constant value by replac ing the data fields lt function gt lt parameters gt by a single number b For a two phase system consisting of a gas g phase and a liquid phase the primary phase will be gas and the capillary pressure parameter will be of the form pe liquid lt function gt lt parameters gt pc which will be used to calculate pe pg pe as a function of aqueous phase saturation c For a two phase system consisting of a NAPL n phase and a aque ous phase the primary phase will be NAPL and capillary pressure parameter will again be of the form pe liquid lt function gt lt parameters gt pc which will then be used to calculate pe pn pe as a function of aqueous phase saturation d For a three phase system consising of gaseous g NAPL n and aqueous L phases the capillary pressure parameter will be of the form for pe gas lt function gt lt parameters gt aqueous lt function gt lt parameters gt pc The parameters for gas defines Pen Py Pn as a function of gas saturation The parameters for aqueous defines pee Pn pe as a function of aqueous phase saturation 4 USING THE COMPUTER PRO
55. isothermal the temperature e uniform and constant temperature XXX APPENDIX b The order of components and phases are important It affects which primary variables will be used depending on which phases are present in a given element See the section on the setting of initial conditions where the state command is described c The primary phase set by primary phase is the phase in which the pressure primary variable P is defined The pressure of the other phases is calculated with respect to this pressure using the capillary pressure relationships See the documentation of material types in rocktab for how these relationships are defined d The wetting phase set by wetting phase is the one which has parti tioning of components between it and the solid phase e For three phase problems the names of phases must be gas NAPL and aqueous with the primary phase set to gas The name liquid can be used instead of aqueous 4 6 Specifying the grid and material property zones The input parameters for the grid and material property zones are genmsh Begin mesh specification coord lt coord type gt Coordinate type down lt down_x gt lt down_y gt lt down_z gt Vector pointing down dx lt dx_1 gt lt dx_2 gt lt dx_nx gt 33 Subdivisions in 33 x direction dy lt dy_1 gt lt dy_2 gt lt dy_ny gt Subdivisions in 33 y direction dz lt dz_1 gt lt dz_2 gt lt dz_nz gt 33 Subdivisions in
56. lvii runiini res instead of file ext res e One can add the statement dstep lt num gt to any of the options ex cept for forcetimes in order to request that the output will be at every lt num gt th time step The outtimes statement is optional in this case If present the counter that keeps track of the dstep option will e reset to zero after every time that hits an output time given in the outtimes statement f Placing the optional statement cumulative in bcflux history will result in the cumulative flux being outputted instead of the in stantaneous flux g One can use the statement index con lt i0 gt lt j0 gt lt KO gt lt i1 gt lt j1 gt lt k1 gt to specify that the connection connects the elements with indices lt i0 gt lt j0 gt lt kO gt and lt i1 gt lt j1 gt lt k1 gt instead of using the connection statement The following is a partial list of output variables Element Variables H Pressure head m for ucsat and ustp only piezometric head if output by piezohead on is present P Primary phase pressure Pa pressure head m for ucsat usip PiezoH Piezometric head m for ucsat Piezo Piezometric head m for usip P lt phase gt Absolute phase pressure Pa Pc lt phase gt Capillary pressure Pa kr lt phase gt Relative permeability S lt phase gt Saturation X lt comp gt lt phase gt C lt comp gt lt phase gt porosity Kx Ky or Kz log10Kx log10Ky
57. mbinative Behie and Forsyth 1984 e Incomplete block LU decomposition with the same ordering options as for the Gauss elimination method For small two dimensional problems the standard Gauss elimination method with natural ordering is best in terms of memory and computa tional costs For larger two dimensional problems the Gauss elimination method with D4 ordering is usually the most efficient For large three dimensional and very large two dimensional problems the preconditioned conjugate gradient method is the most efficient The block Gauss Seidel preconditioning method may be used for problems that are linear such as saturated zone problems For nonlinear problems e g multiphase ones or for problems with large contrasts in permeability a good preconditioning method is the incomplete block LU decomposition with D4 ordering XX APPENDIX 3 7 Reduction of Numerical Dispersion The ability of a numerical method to resolve sharp concentration fronts is reduced by numerical dispersion see Bear and Veruijt 1987 p 323 whenever the grid Peclet number V L D is sufficiently large Numerical dispersion occurs because the discretized equations introduce an artificial diffusive or dispersive flux that is proportional to the grid Peclet number NUFT has two methods for reducing numerical dispersion appearing in the numerical models of the transport equations e Introduction of artificial diffusion which us
58. ng elements Considering the usnt module we desire that the xlii APPENDIX aqueous phase pressure at the centroids of the pond elements be Pe Pg peg He 4 1 where p is atmospheric pressure Therefore the fictitious capillary pressure function for the pond is Pe Pg Pe pegH Se 0 5 peg AZ pond 4 2 Similarly for usip we have De S 0 5 AzZ pona for usip 4 3 Note that unlike a true capillary pressure it can be negative In set ting up the model the linear interpolatory table option for specifying capillary pressure is used to implement the above function For usnt a row of auxilliary elements above the row of ponding ele ments is needed to maintain atmospheric gas pressure Pg Patm in the elements through the use of bctab For ustp a row of completely unsaturated elements S1 O is needed because an element which in this case is a pond element can not become unsaturated unless there is a neighboring element that is unsaturated These added elements are of the same type as those used for the atmospheric boundary elements Surface precipitation with varying head This type of boundary condition is implemented in the same manner as the previous one A row of elements above the ground surface has flux implemented by the srctab option There is a another row above these to maintain a fixed gas pressure Well inflow An element can be used represent a well screen for injection or extr
59. ng subscripts and superscripts A Swp OG VASP WV D Vw D Vw fjs t O80 2 12 where p S and V are transferred from the flow module For the component on the solid we have C p 50 fjs 1 opel 2 13 A source due to decay is expressed by T w T AF NUFT employs the linear isotherm F Kaw 2 4 Transport of multiple components in multiple phases un der isothermal and nonisothermal conditions The module usnt is used for this class of problems It considers the case of NP fluid phases and NC components 2 MATHEMATICAL MODEL vil Since by summing up mass balance equations for all components within a phase we obtain the mass balance equation for that phase NUFT actually solves the NC mass balance equations for the components The mass balance equation for an individual y component in all a phases is of gt Sawa pa R a V 45 Sapa Wl Va D Vol a a 2 dS opal i 1 psF2 2 14 a in which the sum is over all fluid a phases 7 s are sources due to first order growth or negative sink for decay Rj denotes the retardation factor for the y component in the a phase defined as 1 PP Ka RO l da da bS a and the various K4a s may be constant or functions of the temperature NUFT assumes that there is a unique wetting phase amp wet that can form a film on the solid surfaces which can lead to sorption on the surfaces Therefore Kj 0 for all a F awet The phase velo
60. of the pond The concentration of the contaminant is kept at the desired value which can vary in time The concentration of dissolved air in the aqueous phase must also be set The flow of aqueous phase from the pond can be turned on or off at given times by using the phase factor tables The contaminant in the pond can also be shut off using the factor tables or by changing the specified concentration in the pond Ponded contaminant source of varying depth If no water is added to the pond or if the amount of water that is added to the pond changes in time the level of the water in the pond will rise or fall at a rate depending on the infiltration flux which in turn depends on the conditions in the soil The flux of water into the pond is represented using the srctab option in order to specify a time dependent flux into the elements representing the pond The material parameters of the pond are chosen in the same manner as the ponded source of constant depth except that a fictitious capillary pressure function is used to maintain the correct pressure head in the pond as a function of wa ter level The aqueous saturation Sy of the pond elements is used to represent the fraction of each element filled with aqueous liquid as suming vertical equilibrium within each element The water elevation level relative to the height of the centroid of the pond element is given by He Se 0 5 AzZ pond Where AZ pond is the thickness of the pondi
61. oid of elements and variables with double subscripts such as m m are defined at the common boundary of elements m and m The size of a time step is defined as Atn41 tn41 tn The component density p of the mass flux vector perpendicular to the flow area A between elements m and m is expressed in the discretized form pkk 8 OPV intm Tn The relative permeability k m m is evaluated at the upstream value i e Pm Pm Pim m9 2m 2m 3 3 kr Sm if Pm Pm Pint 2m Zm gt 0 kr m m l k Sm otherwise 3 4 The non advective flux dispersive plus diffusive is expressed in the form Dm m OPJ mtm 0p mm Wm Wm 3 5 Lm m 3 NUMERICAL MODEL xi where Lm m is the distance between the centroids of the elements m and m For cylindrical coordinates r y z along the angular y direction the distance Lm m is the length of the arc of constant radius r The flow area in the radial direction between two elements is given by Amtm Tmm AZAY 3 6 Optionally the logarithmic area Tm Tm AzAy Anim In tim Tm 3 7 cal be used which gives the exact steady state solution for radial flow but not for transport The values of oKk H m m correspond to the element that is located upstream relative to the sign of the flow velocity yin The value of Dm m is the harmonic average given by Lint Lm Lit Pm Lm Pm Dn
62. on of the problems presented in Chap although other codes may also be used As stated in the pream ble to that chapter the objectives of presenting these problems and using the NUFT code to solve them are a to gain some experience in the ap Dec 17 95 22 ii APPENDIX plication of numerical models and the use of computer codes b to present solutions to the models discussed in this book and c to investigate how these solutions vary as various input parameters are changed Chap also contains exploratory problems that one may pursue on one s own The NUFT code is available from the International Groundwater Modeling Center IGWMC Golden Colorado USA The version de scribed in this section and distributed with this book runs on IBM com patible personal computers under Microsoft Windows The NUFT code is written in the computer language C It also runs on UNIX workstations The code has been verified by comparison with analytical solutions of flow and transport problems and with numerical solutions by other computer codes Lee et al 1993 This appendix contains some information on the NUFT code and guide lines on how to use it It is written under the assumption that the potential user of NUFT is familiar with the principles of numerical techniques espe cially with those of the method of integrated finite differences It explains how NUFT numerically solves the flow and contaminant transport mod
63. onents 4 USING THE COMPUTER PROGRAM li and phases where the init eqts statement is defined For two phase problems if if the primary phase is not present in an element i e saturation of the primary phase is zero then P becomes the pressure of the non primary phase For three phase problems gas aqueous NAPL phase the primary phase has to be the gas phase If only the gas phase is not present then P becomes the NAPL pressure If both gas and NAPL phases are not present it becomes the aqueous phase pressure S lt phase gt denotes the saturation where lt phase gt is the name of the phase The saturation of the first phase that occurs in the list of phases in init eqt is not set by the user It will be calculated as one minus the sum of the saturation of the other phases set by the user C lt comp gt denotes the concentration where lt comp gt is the name of the component The default units for concentration are in mole fraction If the statement input mass fraction on is present than the units are in mass fractions and the variable name is instead X lt comp gt The concentration in a particular element refers to the component in the first phase with non zero saturation in that element with respect to the ordering of phases set in init eqts If NC is defined as the total number of mass components and NZP is the number phases with non zero saturation in a given element Only the concentrations of the last NC NZP components in the
64. or an isothermal model using usnt module the parameters are porosity lt value gt 33 Kx lt value gt Ky lt value Kz lt value gt Kd lt comp gt lt value gt Kd 33 pe lt phase gt lt function gt lt parameters gt pc kr lt phase gt lt function gt lt parameters gt kr tort lt phase gt lt function gt lt parameters gt tort gt Same as used in ucsat above Sat permeability m2 in the x y z directions The Kx value is used for all directions unless the anisotropic option is set in genmsh Normalized linear iosotherm same definition as for usic above value needed for every mass component Capillary pressure Pa function 3 function See below for which phases are required Relative permeability function function needed for all phases gt Tortuosity factor unitless function one needed for all phases see definition for XXXVI APPENDIX 3 usic module above dispL lt num gt 3 Same as for usic above dispT lt num gt 3 Same as for usic above If a model using the usnt module is nonisothermal the following addi tional parameters are required in addition to those for the isothermal prob lem KdFactor Temperature dependent factor lt comp gt lt function gt lt parameters gt Cunitless multiplying the KdFactor 3 normalized Kd factor Cp lt value gt Solid phase specific
65. or the primary variables The user also determines the size of each FD cell and the total number of cells within the domain Henceforth we shall use the term element or cell interchangeably NUFT is based on the numerical technique called integrated finite differ ences method e g Bear and Verruijt 1987 p 244 It is also referred to as the finite volume method This method allows FD elements of arbitrary polyhedral shapes as long as the line segment connecting the centroids of any pair of neighboring elements is orthogonal to their shared boundary The integrated finite difference method reduces to the standard finite dif ference method when the elements constitute a standard rectangular mesh Here we shall restrict the discussion to FD elements that are identical to those assigned in the standard method of finite differences for rectangular or cylindrical coordinates In cartesian coordinates the elements are 2 d rectangles or 3 d parallelepiped boxes with sides parallel to the axes In cylindrical coordinates they are sectors of an annular cylinder The discretized form of the balance equation at the mth element is Um On pet op ph Att J Amm Cages OPTIA Un OnE 3 2 mm m ENm where superscripts denote time levels and subscripts denote spatial locations with m Nm denoting the set of elements that are neighbors of the mth one Variables with single subscript are those defined at the centr
66. order ap pearing in init eqts are actually used by the model for that element We call these concentrations the independent concentrations for the element The concentrations for the other components are internally computed by the model using the partitioning coefficients as a func tion of the NC NZP independent concentrations the temperature and pressure assuming thermodynamic equilibrium These other concen trations are called the dependent concentrations The user can set dependent concentrations to a negative value The model will give an error message and stop if a negative value is set for an independent concentration Tis the temperature in degrees centigrade For isothermal problems the initial temperature is not set by state but by the following generic statement generic T lt temp gt generic where lt temp gt is the value of temperature in degrees centigrade that will be uniform over the entire model and constant in time lii APPENDIX All primary variables except those stated otherwise above must be set The order in which the primary variables are set in state is does not matter A primary variable can be specified more than once A specification of an initial condition that applies to a given element will override any previous specification That is if the pressure of an element is set in more than one specifications the last one is what the model takes as the initial condition Examples will be given below The variou
67. ot present for isothermal xxxviii APPENDIX pcTemFac Temperature dependent factor multiplying capillary pressure not present for isothermal lt phase gt phaseprop 4 9 Specifying boundary conditions of the first and third type Boundary conditions of the first and third types are implemented through the bctab data record with form given by bctab Begin specification of boundary conditions lt bc name gt Name of this particular boundary condition can be any name range lt range gt lt range gt Range of elements basephase lt phase gt 3 Not present for ucsat 3 or usip modules tables Specify time history of primary variables lt primary var gt lt t0 gt lt valueO gt lt t1 gt lt valuel gt lt primary var gt lt t0 gt lt valueO gt lt t1 gt lt valuel gt tables End tables factor Optional time table of factors which multiplies all component fluxes across this boundary condition lt comp name gt lt tO gt lt v0 gt lt t1 gt lt v1 gt Factor table for component lt comp name gt lt tO gt lt vO gt lt t1 gt lt v1 gt factor phasefactor Optional time table of factors which multiplies all phase fluxes across this boundary condition lt phase name gt lt t0 gt lt v0 gt lt t1 gt lt v1 gt Factor table for phase lt phase name gt lt t0 gt lt v0 gt lt
68. otherwise noted as part of the data input format All characters after a semicolon to the end of a line are comments they are ignored by the program The use of two consecutive semicolons is recommended in order to make comments stand out The MKS system is used throughout Thus for example the following units are used meters for distances Pascal for pressure meter squared for permeability meter per second for hydraulic conductivity and kilo gram per cubic meter for mass density An exception to the above rule is the units for time The units of numeric values that designate time may be optionally represented by appending one of the following characters to the end of the number with no intervening spaces s for seconds m for minutes h for hours d for days M for months and y for years For example 4 2h represents 4 2 hours The default unit no appended character is seconds Temperatures are in centigrade unless otherwise noted A pattern string is a special type of a string of characters which can have unix type wild characters or so that a pattern string can represent an entire class of strings that matches the string pattern The character in a pattern matches any sequence of characters Hence the pattern matches all strings The character in a pattern matches any single character Hence the pattern 7 matches all strings with exactly one character Other Examples
69. ple of primary variables for an air water contaminant system Consider a system with air a water w and a single volatile contaminant c as the components Possible phases are the aqueous phase which is mixture of liquid water and dissolved contaminant and air gaseous phase g which is a mixture of air water vapor and contaminant vapor and NAPL phase n which is a mixture of liquid contaminant with dissolved air and water The number of primary variables is NC 1 4 In the follow ing examples of primary variables the a 3 6 denote distinct phases with a b L g n e If Sa 1 Sg 0 and S 0 NZP 1 then we use pa wi w and T e If Sa Sg gt 0 and Ss 0 then we have pa Sa ws and T e If Sa S8 S gt 0 we have pa Sa Sg and T 3 3 Numerical representation of dispersive flux Consider the dispersive flux between two elements 1 and 2 Fig 3 1 with concentrations Cy and Ch at their centroids Concentrations Ca and Cg are calculated at the vertices A and B by interpolation using centroid values of adjacent elements By applying the finite difference technique to the dispersive flux qf D VC 3 11 and taking the scalar product with respect to n we obtain iy n Dy2 n C2 C L La n Dy2 t Ca CB LAB 3 12 3 NUMERICAL MODEL XV where qf denotes the dispersive flux from element 1 to element 2 and D12 is the dispersion tensor calculated from the flow velocity at the
70. ponent Note that the elements in the specified range can be but do not have to be part of the bc 1st type or bc 3rd type commands set in genmsh in order to set the flow lengths from the elements to the rest of model to zero or unity respectively In this way if one considers the elements as auxiliary elements one can place the boundary condition right at the boundary of the domain consisting of non auxiliary elements The following boundary conditions require special treatment by using auxilliary elements that have special material types Most of them are present in unsaturated problems a Atmosphere In modeling the vadose zone with a volatile component with the usnt model it is usually necessary to have elements which rep resents the atmosphere above the ground surface in order to model the gas phase transport between the vadose zone and atmosphere This boundary condition is modeled by a row of auxiliary elements at the top of the model implemented as a boundary condition of third type with specified gas pressure concentration of water vapor and if the model is nonisothermal temperature The gas saturation of the ele ments is set to be one Since the resistance to gas flow is negligible compared to the soil a permeability that is two orders of magnitude greater than the maximum permeability of the soil surface can be used The half thickness of the elements is set to the estimated thickness of the mixing boundary layer over the ground
71. raction at specified points in the aquifer pumping and recharging wells For essentially horizontal flow in a phreatic aquifer the variable is the water table elevation h h x y t and the flux and mass balance equations are Oh Sua V Q Qulz y t R x y t Q h n K Vh 2 7 2 MATHEMATICAL MODEL v where 7 7 a y denotes the elevation of the aquifer s bottom and R represents a distributed source e g due to natural replenishment from pre cipitation In both confined and phreatic aquifer balance equations we may add on the r h s source terms that represent leakage through aquitards into or out of the aquifers Such a term will take the form hem h w Cr where q denotes leakage into the aquifer as volume per unit area per unit time Res denotes the head in the underlying or overlying aquifer and c B K denotes the resistance of the aquitard of thickness B and hydraulic conductivity K Phreatic surface as a boundary In Subs we have discussed the conditions along a phreatic surface which serves as the upper boundary of a three dimensional porous medium domain It was shown there that for the case of a phreatic surface with accretion the condition of equality of fluxes normal to this boundary takes the form of 72 Oh Pett KVh N V h 2 2 8 where we have omitted the subscript w N NVz denotes the rate of accretion and denotes the effective porosity
72. rdt the time step size is cutback and the time step is redone The size of the time step for the next step is based on how close the actual time step was to the tolerance Default off rmstolerdt lt num gt Tolerance for the root mean square of the estimate of the relative time truncation error Default 0 1 rms NR conv lt on or off gt Turns on or off Newton Raphson convergence based on root mean square norm of the residual vector R Default off rmstolerconv lt num gt Tolerance for Newton Raphson convergence based on root mean square norm of the residual vector R Default 0 001 5 Running the Exercises in the Book 5 1 Problem 9 1 1 5 2 Problem 9 2 1 Aziz K and A Settari Petroleum Reservoir Simulation Applied Science Publishers London 1979 Bear J and J J Nitao On equilibrium and the number of degrees of freedom in modeling transport in porous media Transport in Porous Media vol 18 no 2 151 184 1995 Bear J and A Verruijt Modeling Groundwater Flow and Pollution D Reidel Pub Co Dordrecht The Netherlands 1987 65 RUNNING THE EXERCISES IN THE BOOK lxii Behie G A and P A Forsyth Jr Incomplete factorization methods for fully implicit simulation of enhanced oil recovery SIAM J Sci Stat Computing 5 543 561 1984 Brooks A and T J R Hughes Streamline upwind petrov galerkin formu lations for convection dominated flows with particular emphasis on the incompressible navier
73. repared the pro gram is run by typing nuft filename inp in a DOS prompt window where filename inp stands for the name of the input file of the considered problem While running the program will output to the window the following in formation at each time step is time step number dt length of previous time step ndt length of current time step iter number of Newton Raphson iterations the option to output the primary variable that is restricting the time step In addition the program has size For example for the us1p module the output line dt restricting Sl 0 35 dSl 0 15 H 1 1 2 indicates that the time step was restricted by a change in liquid saturation S1 from the previous time step by 0 15 to the new value of 0 35 at element Ri 1 1 2 Note that NUFT refers to elements by their names which take the form lt prefix gt lt i gt lt j gt lt k gt where lt prefix gt is usually the name of a region here R1 and lt i gt lt j gt lt k gt here 1 1 2 are the i j k indices of the element within the grid As requested in the input file NUFT prepares an output file using the above time step information The name of the output file has the form lt run gt out where lt run gt is the the run name which is identical to the file name of the input file For example the run with the name p911 inp will have an output file called p911 out NUFT will also create other output files upon request In parti
74. rgy flux W m in the phase is output Specific conductive thermal flux W m Total flux of phase kg s Total flux of a component kg s in all phases if lt comp gt is energy then total energy flux W in all phases including solid phase is output Total flux of a component kg s in a phase if lt comp gt is energy then total energy flux W in the phase is output Total conductive thermal flux W The following parameters control various output options Coutput prefix lt prefix gt If present all output files will consist of this prefix together with the suffix of the difference files For example set ting Coutput prefix runi will mean that the output files will be runt out for the main output file and runi nvc and runt nvh 4 USING THE COMPUTER PROGRAM lix for the nview files If not present the model will use the first part of the name the input file until the first period if any in its name as the output prefix Default not present print mat bal lt on or off gt Turns on or off the printing of material balances Valid options are on or off Default off print pcg lt on or off gt Turns on or off the printing of preconditioned conjugate gradient information Valid options are on or off Default off tty output step lt num gt Output time step information every this many time steps Default 1 which means at every time step 4 15 Numerical Method Options The following
75. rsion on gt APPENDIX is employed Optionally read in initial conditions from another run not used if state is employed Specify output variables and times Specify component properties present only for usic and usnt modules Specify fluid phase properties present only for usnt module Present only for isothermal problem using usnt module used for setting initial temperature to lt temp gt in degrees centigrade If present units of concentrations in input file is in mass fractions If not present concentrations are in mole fractions Applies only to ucsat and usip modules if present the initial conditions for the variable H are specified as piezometric head instead of pressure head If present the model implements hydrodynamic dispersion use only for usic and usnt modules numerical method options lt module name gt gt End module input parameters To use flow and transport modules together such as ucsat for flow and usic for transport two sets of statements of the kind presented above would be required each containing input parameters that are specific to one of the models In addition there exists a common set of data which should also be included This is implemented by the command common Place here common input parameters such as common 3 genmsh data The following item is included in a flow model 4 USING T
76. s by using the overwrite restart option which has exactly the same syntax and options as the state option described earlier 4 14 Specifying output options output contour Output element variables 3 to nvc NVIEW input file for contouring and vectors variables lt vari gt lt var2 gt 3 Names of element variables outtimes lt ti gt lt t2 gt List of output times contour elem history Output element variables 3 to nvh NVIEW input file for time history plots variables lt vari gt lt var2 gt 3 Names of element variables range lt range gt lt range gt Range of elements outtimes lt ti gt lt t2 gt List of output times elem history connect history Output connection variables 3 to nvh NVIEW input file for time history plots variables lt vari gt lt var2 gt 3 Names of connection variables connection lt elemO gt lt elemi gt Names of elements of the 33 connection outtimes lt ti gt lt t2 gt List of output times connect history lvi APPENDIX bcflux history Output boundary condition fluxes to nvh file variables lt var1 gt lt var1 gt Names of fluxes name lt bc name gt Name of boundary condition specified in bctab outtimes lt t1 gt lt t2 gt List of output times bcflux history restart Output primary variables to a file for restart runs Couttimes lt t1 gt
77. s initialization methods are as follows a lt primary variable gt by key lt range gt lt value gt lt range gt lt value gt The primary variable is set to the lt value gt for all elements that are specified by the associated element range string pattern lt range gt If an element falls within more than one string pattern the value asso ciated with the last pattern that matches that element is used For example P by key 1e5 f 2e5 first specifies that all elements have pressure set to 1e5 and then over writes the elements that start with the character f to the value 2e5 b lt primary variable gt by xtable range lt elem range gt lt x0 gt lt value0 gt lt x1 gt lt valuel gt lt primary variable gt by ytable range lt elem range gt lt yO gt lt valueO gt lt y1 gt lt value1 gt lt primary variable gt by ztable range lt elem range gt lt z0 gt lt value0 gt lt zi gt lt valuel gt Any one or more of the above methods can be used The methods lt by xtable gt lt by ytable gt and lt by ztable gt sets the initial condi tions from linear interpolation of a table with respect to the coordi nates i yj Zk of the centroid of each element respectively The centroid of an element with indices i 7 k is given by x Si Azi y Dia Ay Z Si Azy where Azry Ayy Az are the grid subdivisions of the elements as set in the genmsh sta
78. sources distributed and point sources their loca tions and strength Here we include the locations and rates of pumping and 4 USING THE COMPUTER PROGRAM xxi recharging wells We also identify the initial conditions with respect to the relevant vari ables Zones of common initial conditions may be specified Sometimes we may wish to start from an initial steady state the details of which are a priori unknown We obtain this steady state regime by running an initial ization problem in which we fix boundary conditions and start from any arbitrary initial conditions After a sufficiently long time depending on how close are the suggested initial conditions to the ultimate steady state the system will reach steady state This steady state regime is then used as in put to the considered problem NUFT has the option of defining the output from the initialization problem as a lt filename gt rst standing for restart and inserting this lt filename gt rst file as input for initial conditions in the input file filename inp We continue by identifying all the relevant physical processes actually also the chemical and biological ones but here we focus on the physical processes only Some important questions to be asked are e Does the considered problem refer to a saturated domain If not which fluid phases other than water or an aqueous fluid are involved in the problem e Which if any of the components are volatile e
79. stribution of these state variables at specified times NUFT consists of four different models within a single code a ucsat unconfined aquifer saturated flow model b usip unsaturated single phase flow model c us1c unsaturated single component transport model d usnt fully coupled unsaturated multiple phases multiple compo nents model with isothermal and nonisothermal options LON INN All models share a common set of internal utility routines e g for input output and numerical algorithms For a problem of flow in a saturated zone or a problem of single phase flow in the unsaturated zone rather than using the usnt model it is compu tationally more efficient to use the ucsat and usip models respectively and then to use the usic model for solving the contaminant transport problem running alongside ucsat or us1p The usnt model has the option to solve the energy balance equation coupled to the mass balance equations thus facilitating the solution of non isothermal flow and transport problems Of particular importance espe cially for nonisothermal problems is NUFT s ability to handle properly the appearance or disappearance of any phase due to condensation and evapo ration 2 Mathematical Model Actually all the mathematical models especially the p d e s solved by NUFT have already been presented in the appropriate chapters of the book To facilitate the use of NUFT let us rewrite them in the forms that
80. t 3 8 Note that the right hand side of 3 10 is evaluated at the n 1 th time step which means that a set of equations usually nonlinear must be solved at each time step This method of time discretization is called implicit in time Evaluation at the n th time step is called explicit in time NUFT has the option to use either method The implicit in time method is more costly to compute at each time step but often much larger time steps can be taken because of its increased numerical stability If the r h s of 3 2 defined at time levels n and n 1 is denoted by F and F respectively NUFT also allows a linear weighting for the r h s in 3 2 af 1 a F where a 0 lt a gt 1 is a weighting constant Note that a 1 is fully implicit in time and a 0 is explicit in time The option a 0 5 is the so called Crank Nicholson method which gives second order error in time discretization whereas the other two methods are only first order in time However this method can lead to oscillations in time The option a 0 6 still gives good accuracy but is more stable and is recommended especially for transport problems If we write 3 2 in the residual form R v 0 3 9 xii APPENDIX where v denotes the solution at time level n 1 then the Newton Raphson method Aziz and Settari 1979 applied to this set of nonlinear equations is JO vetD yh Rv v 0 1 3 10 where the superscript
81. t dir gt The flux at element m equals the specified flux multiplied by the factor Am X m Am Where m runs over all elements specified in the range of the particular source term and Am are element flow areas Since the flow area of element is different depending on its sides the direction of the flow area of an element is specified by the lt dir gt parameter Valid options for lt dir gt are x y or z corresponding to the z y and z directions respectively In cylindrical coordinates these directions correspond to the r and z directions and the flow area in r is set to the average of the two flow areas in the r direction except for the elements with i l and nz in which case only a single flow area is used that of the flow area common to the neighboring element in the r direction Multiply flux by flow area of the particular element mult by area lt dir gt The flux at element m equals the specified flux multiplied by the flow area A of the element in the lt dir gt direction where lt dir gt is the same as the allocate by area option above and the same conventions for determining flow area apply Multiply flux by factor which depends on a group of elements allocate by element lt range gt lt fac gt lt range gt lt fac gt allocate by element If an element is not in the specified ranges the flux value is that computed from the table Otherwise if an element is within a range the fl
82. tement described earlier The coordinate values in the table must be strictly increasing order The minimum coordinate value in the table must be less than or equal to the coordinates of the centroids of all elements that are specified by range and the maximum coordinate must be greater than or equal to the coordinates of the centroids of all elements specified by the range 4 USING THE COMPUTER PROGRAM liii The values obtained from the table are set only for the elements in the range specified by the range statement The range statement is optional if not present all elements in the model will be set As mentioned before the index range statement can be used in place of the range statement to specify the range of elements Sometimes it is useful to generate a table of primary variables from a one dimensional initialization run NUFT can output primary vari ables in table format See the NUFT Reference Manual for details c lt primary variable gt by ijk file lt file name gt The model will read initial conditions from a file The form of the file must be lt var name gt lt 10 gt lt j0 gt lt k0 gt lt value0 gt lt i1 gt lt j1 gt lt k1 gt lt value1 gt end where lt var name gt is the character followed immediately by the name of the primary variable Each subsequent lines of the file contains the 2 j k index of an element and its initial condition Every active i e non null element must have a line Any
83. the 3 NUMERICAL MODEL xix auxilliary element which multiplies the Fickian diffusion coefficient and the dispersivities If desired setting the tortuosity factor and dispersivity of the auxiliary element to zero will eliminate the dispersive flux in the above boundary condition Some types of boundary conditions especially for unsaturated problems require specialized techniques in their implementation and will be described later in the section on the use of the computer model 3 6 Solving the system of linear equations Equation 3 10 requires the solution of a system of linear equations at each iteration NUFT can implement several different methods for solving the system of linear equations Which method is the most efficient in terms of computational memory and CPU time depends on the number of elements and on the dimensionality of the problem The methods included in NUFT are e Gauss elimination with various orderings of the components of the so lution vector including natural ordering D4 ordering Aziz and Settari 1979 and ordering by the reverse Cuthill McKee bandwidth minimiza tion algorithm Cuthill and McKee 1969 Preconditioned conjugate gradient method using the ORTHOMIN algo rithm Behie and Forsyth 1984 which is an iterative method This method requires that the matrix be preconditioned by multiplying it by an approximate inverse matrix The available preconditioning options are e Block gauss seidel co
84. tion e An unsaturated element is switched to saturated in the course of a simu lation if the element s saturation becomes greater than or equal to unity e A saturated element is switched to unsaturated if the pressure head in that element becomes negative An important exception in both ucsat and us1p modules is that an ele ment is not switched to unsaturated conditions if there exists no neighboring element that is unsaturated For the usnt module let NC and NP denote the number of components and number of phases respectively It can be shown ref Bear and Nitao 3 NUMERICAL MODEL xiii 1995 that there are NC primary variables if the problem is isothermal and NC 1 primary variables if the problem is nonisothermal To determine the primary variables let NZP denote the number of non zero phase saturations the actual number of phases that exist in an element i e for the saturation Sa of an a phase we have Sa gt 0 a 1 NZP Let NC denote the number of components and let a denote some phase such that Sy gt 0 Then the following are used e pressure P of one selected phase e NZP 1 saturations Sa amp 1 NZP 1 e NC NZP mole fractions n y 1 NC NZP and e temperature T if the case is nonisothermal At t 0 we know the number NZP and the primary variables are selected For t gt 0 the code will determine the primary variables as the number NZP varies always using the last N
85. ux at that element equals the value computed from the table multiplied by the respective factor lt fac gt Note that an 4 USING THE COMPUTER PROGRAM xlvii element can be in more than one range in which the flux will be multiplied by more than one factor 4 11 Specifying wells Although it is possible to use srctab to implement wells it is often more convenient to use the special well option It is of the form setwells lt well name gt Name of well used for output purposes rw lt num gt Radius of well bore m re lt num gt Radius of pressure support m lt well type gt lt well type gt is one of following 33 pres inj specified pressure injector 33 pres prod specified pressure producer 33 producer can only have flow out of formation injector only flow into 33 the formation see below for rest of input data lt well type gt intervals Specifiy well intervals intervals must be in order of topmost first i e highest actual elevation first setint Specify parameters for a single interval index lt i gt lt j gt lt k gt Specify i j k index of element zcoord lt num gt Optional specifies z coord of 3 midpoint of interval 3 w r t mesh coord system 3 default is z coord of element center length lt num gt 3 Optional specifies length of screened portion of interval default is z thickness of the completed 3 element
86. y a If the pressure head H is set to a value greater than zero then then the model uses pressure head H as the primary variable The value of S1 is set to the maximum saturation Smax given in rocktab b If S1 is set to a value greater than or equal to Smax then the model sets pressure head H to be zero and uses pressure head H as the primary variable The variable S1 is set to be exactly equal to Smax c If pressure head H is set to a value less than zero then the model sets S1 to Smax and pressure head H to zero and uses S1 as the primary variable d If none of the above applies that is if S1 is less than Smax and pressure head H is set to zero then the model uses S1 as the primary variable Note that in the above H is pressure head If piezometric head is specified the model converts to pressure head before applying the above rules For the ucsat module the following rule applies If S1 is less than unity then pressure head H is set to zero and the primary variable is S1 Otherwise S1 is set to unity and the primary variable is pressure head H The primary variable for the usic module is C lt comp gt which denotes the concentration in any units of the transported component called lt comp gt whose name is set by init eqts For the usnt module the following primary variables must be set a P denotes the pressure of the primary phase in Pascals The defini tion of primary phase is given in the section on names of comp
87. y difference between nuftini and nuft The program nuftini is used to make an initialization run to obtain steady state conditions Thus one can have a single input file with input used only during initialization inside the ifdef INIT statements and non initialization input inside the ifndef INIT statements 4 4 General Content of an Input File Following is the general content of an input file lt module name gt Set lt module name gt to be name of module usip ucsat etc init eqts Specify the names of components wae and phases isothermal or not init eqts 3 not present for ucsat usip usic modelname lt name gt 3 Name of this model can be 3 any name necessary because NUFT has option to run multiple modules time lt t0 gt 3 Initial time dt lt dt0O gt Initial time step stepmax lt N gt Maximum number of time steps genmsh genmsh Specify grid and material zones rocktab rocktab Specify material properties bctab bctab Specify 1st and 3rd type boundary conditions srctab srctab Specify 2nd type boundary 33 conditions and sources state state Specify initial conditions 3 not used if read restart xxviii read restart read restart output output compprop compprop phaseprop phaseprop generic T lt temp gt generic input mass fraction on input by piezohead on dispe

APPENDIX - Index of

Contents

Download Pdf Manuals

Related Search

Related Contents