PARAM

PARAM introduces computation parameters, time stepping information, and default initial conditions.

Record PARAM.1

Free format for 34 parameters, or Format (2I2, 3I4, 24I1, 3E10.4, 2I5).

NOITE, KDATA, MCYC, MSEC, MCYPR, (MOP(I), I = 1, 24), Diffu0, TEXP, BE,

SAVE_Fq, TimeS_Fq

NOITE specifies the maximum number of Newtonian iterations per time step (default is 8)

KDATA specifies amount of printout (default is 1).

0 or 1: print a selection of element variables.

2: in addition, print a selection of connection variables.

3: in addition, print a selection of generation variables.

If the above values for KDATA are increased by 10, printout will occur after each Newton-Raphson iteration (not just after convergence). If negative values are used, printout for the selection of variables will be generated in the file format of TECPLOT. The amount of printout is given by the absolute value of KDATA.

MCYC maximum number of time steps to be calculated.

MSEC maximum duration, in CPU seconds, of the simulation (default is infinite).

MCYPR printout will occur for every multiple of MCYPR steps (default is 10000).

MOP(I) I = 1,24 allows choice of various options, which are also documented in printed output from a TOUGH4 run.

MOP(1) if unequal 0, a short printout for non-convergent iterations will be generated.

MOP(2) through MOP(6) generate additional printout in various sub-routines, if set unequal 0. This feature should not be needed in normal applications, but it will be convenient when a user suspects a bug and wishes to examine the inner workings of the code. The amount of printout increases with MOP(I) (consult source code listings for details).

MOP(2) TimeStepping (main subroutine).

MOP(3) In subroutines for mass accumulation and flux calculations.

MOP(4) No longer used in TOUGH4.

MOP(5) Subroutines for equation of state calculations.

MOP(6) No longer used.

MOP(7) if unequal 0, a printout of input data will be provided.

MOP(8) if ISOT < 0, chooses the option for reducing fracture-matrix interface area (See Appendix D).

MOP(9) determines the composition of produced fluid with the MASS option (see GENER, below). The relative amounts of phases are determined as follows:

0: according to relative mobilities in the source element.

1: produced source fluid has the same phase composition as the producing element.

MOP(10) chooses the interpolation formula for heat conductivity as a function of liquid saturation (Sl)

0: C(Sl) = CDRY + SQRT(Sl)* [CWET - CDRY]

1: C(Sl) = CDRY + Sl * (CWET - CDRY)

MOP(11) determines evaluation of mobility and permeability at interfaces.

0: mobilities are upstream weighted with WUP (see PARAM.3), permeability is upstream weighted.

1: mobilities are averaged between adjacent elements, permeability is upstream weighted.

2: mobilities are upstream weighted, permeability is harmonic weighted.

3: mobilities are averaged between adjacent elements, permeability is harmonic weighted.

4: mobility and permeability are both harmonic weighted.

MOP(12) determines interpolation procedure for time dependent sink/source data (flow rates and enthalpies).

0: triple linear interpolation; tabular data are used to obtain interpolated rates and enthalpies for the beginning and end of the time step; the average of these values is then used.

1: step function option; rates and enthalpies are taken as averages of the table values corresponding to the beginning and end of the time step.

2: rigorous step rate capability for time dependent generation data.

A set of times $t_i$ and generation rates $q_i$ provided in data block GENER is interpreted to mean that sink/source rates are piecewise constant and change in discontinuous fashion at table points. Specifically, generation is assumed to occur at constant rate $q_i$during the time interval [ $t_i$ , $t_{i+1}$ ], and changes to $q_{i+1}$ at $t_{i+1}$. Actual rate used during a time step that ends at time t, with $t_i$≤ t ≤ $t_{i+1}$, is automatically adjusted in such a way that total cumulative exchanged mass at time t is

$Q(t) = \int _0^tqdt'=\sum_{j=1}^{i-1}q_j(t_{j+1}-t_j)+q_i(t-t_i)$ (8-4)

is rigorously conserved. If also tabular data for enthalpies are given, an analogous adjustment is made for fluid enthalpy, to preserve $\int qhdt$.

MOP(13) determines content of INCON and SAVE files.

0: standard content.

1: writes user-specified initial conditions to file SAVE, including permeabilities.

2: reads parameters of hysteresis model from file INCON.

MOP(14) (void)

MOP(15) determines conductive heat exchange with impermeable geologic formations using semi-analytical methods (see Section "Semi-Analytical Conductive Heat Exchange").

0: heat exchange is off.

1: linear heat exchange between a reservoir and confining beds is on (for grid blocks that have a non-zero heat transfer area; see data block ELEME). Initial temperature in cap- or base-rock is assumed uniform and taken as the temperature with which the last element in data block ELEME is initialized. If a rock named "QLOSS" exists, the initial temperature can also be inputted through POR in record ROCKS.1 for this rock.

2: linear heat exchange between a reservoir and confining beds is on. Initial temperature for the confining bed adjacent to an element that has a non-zero heat transfer area is taken as the temperature of that element in data block INCON.

Heat capacity and conductivity of the confining beds are specified in block ROCKS for the particular domain with a rock name "QLOSS" (if it exists in ROCKS data section) or to which the very last volume element in data block ELEME belongs. Thus, if a semi-analytical heat exchange calculation is desired, the user would append an additional dummy element of a very large volume at the end of the ELEME block, and provide the desired parameters as initial conditions and domain data, respectively, for this element. It is necessary to specify which elements have an interface area with the confining beds, and to give the magnitude of this interface area. This information is input as parameter AHTX of volume element data in block ELEME (data lot 11). Volume elements for which a zero-interface area is specified will not be subject to heat exchange.

5: radial heat exchange between fluids in a wellbore and the surrounding formation is on. Constant well and formation properties are given in a material named QLOSS with the following parameters:

DROK: Rock grain density [kg/m3] of formation near well

POR: Well radius [m]

PER(1): Reference elevation [m]; specify Z coordinate in block ELEME, Columns 71–80

PER(2): Reference temperature [°C]

PER(3): Geothermal gradient [°C/m]

CWET: Heat conductivity [W/kg °C] of formation near well

SPHT: Rock grain specific heat [J/kg °C] of formation near well

6: radial heat exchange between fluids in a wellbore and the surrounding formation is on. Depth-dependent well and formation properties (depth, radius, temperature, conductivity, density, capacity) are provided on an external file named radqloss.dat with the information in the following format:

On the first line, NMATQLOSS is the number of elevations with geometric and thermal data, and NMATQLOSS lines are provided with the following data in free format: elevation [m], well radius [m], initial temperature [°C], CWET, DROK, SPHT. Between elevations, properties are calculated using linear interpolation.

RESTART is possible for linear heat exchange between a reservoir and confining beds (MOP(15)= 1 or 2), but not for radial heat exchange (MOP(15)= 5 or 6). The data necessary for continuing the linear heat exchange calculation in a TOUGH4 continuation run are written onto a text file named TABLE. When restarting a simulation run, this file has to be present in the working folder. If file TABLE is absent, TOUGH4 assumes that no prior heat exchange with confining beds has taken place.

MOP(16) provides automatic time step control.

0: automatic time stepping based on maximum change in saturation.

1: automatic time stepping based on number of iterations needed for convergence.

>1: time step size will be doubled if convergence occurs within ITER ≤ MOP(16) Newton-Raphson iterations. It is recommended to set MOP(16) in the range of 2 - 4.

MOP(17) handles time stepping after linear equation solver failure.

0: reduce time step after linear equation solver failure.

>0: no time step reduction despite linear equation solver failure.

MOP(18) selects handling of interface density.

0: perform upstream weighting for interface density.

>0: average interface density between the two grid blocks. However, when one of the two phase saturations is zero, upstream weighting will be performed.

MOP(19) switch used by different EOS modules for conversion of primary variables. It is effect only for EOS3, EOS4, and EOS6. See Section "Inputs for Initial Conditions" for details.

MOP(20) switch for vapor pressure lowering in EOS4; MOP(20)=1 optionally suppresses vapor pressure lowering effects. Users may use IE(24) to turn on the vapor pressure lowering for other EOS modules.

MOP(21) (void)

MOP(22) (void)

MOP(23) (void)

MOP(24) determines handling of multiphase diffusive fluxes at interfaces.

0: harmonic weighting of fully coupled effective multiphase diffusivity.

1: separate harmonic weighting of gas and liquid phase diffusivities.

DIFFU0 base diffusion coefficient in gas phase (not used in current version).

TEXP parameter for temperature dependence of gas phase diffusion coefficient.

BE (optional) parameter for effective strength of enhanced vapor diffusion; if set to a non-zero value, will replace the parameter group for vapor diffusion (see Eq. (4-9)).

SAVE_Fq The frequency for writing out SAVE file. SAVE file will be written every SAVE_Fq time step.

TimeS_Fq The frequency for writing time series output. Time series output file will be written every TimeS_Fq time step, Default is 1.

Record PARAM.2

Free format for 8 parameters, Format (4E10.4, A5, 5X, 3E10.4)

TSTART, TIMAX, DELTEN, DELTMX, ELST, GF, REDLT, SCALE

TSTART starting time of simulation in seconds (default is 0).

TIMAX time in seconds at which simulation should stop (default is infinite).

DELTEN length of time steps in seconds. If DELTEN is a negative integer, DELTEN = -NDLT, the program will proceed to read NDLT records with time step information. Note that - NDLT must be provided as a floating point number, with decimal point.

DELTMX upper limit for time step size in seconds (default is infinite)

ELST The name of an element to be tracked during time-stepping and its key thermophysical parameters will be written out in the log file. For the formatted input, only an element name with 5 characters is allowed.

GF magnitude (m/sec2) of the gravitational acceleration vector. Blank or zero gives "no gravity" calculation.

REDLT factor by which time step is reduced in case of convergence failure or other problems (default is 4).

SCALE scale factor to change the size of the mesh (default = 1.0).

Record PARAM.2.1, 2.2, etc.

Free format for as more as 20 parameters, or Format (20E10.4)

(DLT(I), I = 1, -NDLT)

DLT(I) Length (in seconds) of time step I.

This set of records is optional for DELTEN = - NDLT, a negative integer. Input up to 100 time-step data. Multiple-line input is allowed, but each line inputs at least 1 data and at most 20 data. If the number of simulated time steps exceeds the number of DLT(I), the simulation will continue with time steps determined by automatic time step control (MOP(16)).

Record PARAM.3

Free format for 11 parameters, Format (11E10.4)

RE1, RE2, U, WUP, WNR, DFAC, FOR, AMRES, prvCR1, prvCr2, prvCR3

RE1 convergence criterion for relative error, see Appendix B (default= 10-5).

RE2 convergence criterion for absolute error, see Appendix B (default= 1).

U (void).

WUP upstream weighting factor for mobilities and enthalpies at interfaces (default = 1.0 is recommended). 0 ≤ WUP ≤ 1.

WNR weighting factor for increments in Newton/Raphson - iteration (default = 1.0 is recommended). 0 < WNR ≤ 1.

DFAC increment factor for numerically computing derivatives (default value is DFAC = 10-k/2, where k, evaluated internally, is the number of significant digits of the floating point processor used; for 64-bit arithmetic, DFAC ≈ 10-8).

FOR weighing factor determining the level of implicitness (default =1.0).

AMRES Not be used in TOUGH4

prvCR1 convergence criterion for primary variable change, pressure (default=100.0).

prvCR2 convergence criterion for primary variable change, saturation, or mass fraction (default=0.001).

prvCR3 convergence criterion for primary variable change, temperature (default=0.1).

Record PARAM.4 (optional) introduces the system phase-state index for the default primary variable set inputted through PARAM.5. It is important only when the program cannot figure out its phase state through the primary variables.

Free format for 1 parameter, or Format (1I2)

DefaultStateIndex phase state index, can be found from tables in Process Modeling.

Record PARAM.5 introduces a set of primary variables which are used as default initial conditions for all grid blocks that are not assigned by means of data blocks INDOM or INCON.

Free format for as more as 20 parameters, or Format (4E20.13)

DEP(I), I = 1, NumPriV

The number of primary variables, NumPriv, is determined by the mass and energy component number of the flow system. If TOUGH3 initial condition is used as input, the NumPriv equals NKIN+1. NKIN is given in data block MULTI for special assignments. For formatted input, when more than four primary variables are used, more than one line (record) must be provided. Different sets of primary variables are in use for different EOS modules; see the description of EOS modules.

If the phase state index has not been inputted with record PARAM.4, user can introduce it by inserting -DefaultStateIndex*10 into any location of this record. TOUGH4 will understand that it is for the DefaultStateIndex, if the data is less or equal -10.

Used in: All EOS modules

Example:

PARAM

9, 2, 3000, ,1000, 1, 0*3, 3, 0*8, 2, , 4

// , , MCYC, ,MCYPR, , MOP(2-4), , MOP(6-13), , , MOP(16)

0.0, 27000.0, -1., 2000., , 9.8060 //TSTART, TIMAX, DELTEN, DELTMX, , GF

20. //DLT(1)

1.E-5, 1.E00, , , , 1.e-06, , , 1000, 0.01, 0.1 // RE1, RE2, , , , DFAC, , , prvCR1, prvCr2, prvCR3

2 // phase state index.

1.0133e5, 0.00, 0.0, 22.0 // default initial conditions.

The lase two lines can also be combined into one line:

1.0133e5, 0.00, 0.0, 22.0, -20.0