Initial and Boundary Conditions#
Initial Conditions and Restarting#
Flow systems are initialized by assigning a complete set of primary thermodynamic variables to all grid blocks into which the flow domain is discretized. Various options are available in a hierarchical system, as follows. During the initialization of a TOUGH3 run, all grid blocks are first assigned to default thermodynamic conditions specified in data block PARAM.4. The defaults can be overwritten for selected reservoir domains by assigning domain-specific conditions in data block INDOM. These in turn may be superseded by thermodynamic conditions assigned to individual grid blocks in data block INCON. A disk file INCON written to the same specifications as data block INCON may also be used.
Different EOS modules use different sets of primary variables, and some EOS modules allow different choices of primary variables for initialization. For a given EOS module, the primary variables used depend on the fluid phase composition. During phase change, primary variables will be automatically switched from one set to another. In multiphase flow systems, therefore, different grid blocks will in general have different sets of primary variables, and must be initialized accordingly.
For many applications, special initial conditions are needed, such as gravity-capillary equilibrium, or steady state corresponding to certain mass and heat flows. This can be realized by performing a series of TOUGH3 runs, in which thermodynamic conditions obtained in one run, and written to disk file SAVE, are used as initial conditions in a subsequent continuation run. For example, in a geothermal reservoir simulation, a first run may be made to obtain hydrostatic pressure conditions. These may subsequently be used as boundary conditions in a second run segment to simulate undisturbed natural state conditions with throughflow of mass and heat. This could be followed by a third run segment with fluid production and injection. Restarting of a TOUGH3 run is accomplished by providing file SAVE generated in a previous run as file INCON for initialization. Usually additional (often minor) adjustments will be made for a restart. For example, different specifications for the number of time steps and desired printout times may be made. Some editing of the SAVE and/or MESH files may be needed to implement boundary conditions on certain grid blocks. Then, previously calculated pressures on those grid blocks can serve as boundary conditions. In a continuation run, simulation time and time step counters may be continuously incremented, or they may be reset to zero. For example, the latter option will be used when simulating production and injection operations following preparation of a “natural” initial state, which may correspond to a large simulation time.
As far as the internal workings of the code is concerned, there is no difference between a fresh start of a simulation and a restart.
The only feature that makes a simulation a continuation run is that the INCON data were generated by a previous TOUGH3 run, rather than having them explicitly provided by the user.
File SAVE always ends with a data record with +++ in the first three columns, followed by one record with restart information (time step and iteration counters, simulation time).
To reset all counters to zero, users should simply replace the +++ with a blank record when using SAVE as file INCON for another TOUGH3 run.
Neumann Boundary Conditions#
Neumann conditions prescribe fluxes of mass or heat crossing boundary surfaces. A special case of Neumann boundary conditions is “no flux”, which in the integral finite difference framework is handled simply by not specifying any flow connections across the boundary. More general flux conditions are prescribed by introducing sinks or sources of appropriate strength into the elements adjacent to the boundary. Time-dependent Neumann conditions can be specified by making sink and source rates time dependent.
Dirichlet Boundary Conditions#
Dirichlet conditions prescribe fixed thermodynamic conditions, such as pressure, temperature, or saturation. Dirichlet conditions can be implemented by assigning very large volumes (e.g., \(V\) = 1050 m3) to grid blocks adjacent to the boundary, so that their thermodynamic conditions do not change at all from fluid or heat exchange with finite-size blocks in the flow domain. For those elements with very large volumes, no mass or energy balance equations are set up, and their primary thermodynamic variables are not included in the list of unknowns. In addition, a small value (such as \(D`\) = 10-9 m) should be specified for the nodal distance of such blocks, so that boundary conditions are in fact maintained in close proximity to the surface where they are desired, and not at some distance from it.
TOUGH3 no longer supports the so-called “inactive” elements allowed in TOUGH2. If an element with a zero or negative volume is used to indicate that the element and all subsequent elements are to be used for Dirichlet boundary conditions, then the code will internally assign a very large volume for that and all subsequent elements (to be backward compatible with TOUGH2) and print out a warning message. However, if a dummy element with no volume and no material name is used just to separate boundary elements, the code will generate an error message and stop the simulation. In TOUGH3, the user must specify a material name for every element.
Time-dependent Dirichlet conditions can be implemented using two different approaches. The first method consists of placing appropriate, large sinks or sources in boundary elements with large volumes (Moridis, 1992). As an example, consider a laboratory experiment reported by Kruger and Ramey (1974) that involved flashing (vaporizing) flow from a sandstone core, with a time-dependent gas pressure boundary condition
maintained at the outflow end. According to the ideal gas law, the pressure behavior of Eq. (42) can be associated with a time-dependent gas inventory of mass \(M_b\) in a volume \(V\) as follows:
where \(m\) is the molar weight of the gas, \(R\) the universal gas constant, and \(T\) absolute temperature. The required time dependence of \(M_b\) can be realized by means of a sink/source rate of
Therefore, the desired boundary condition can be implemented by means of a grid block with volume \(V\) that is initialized in single-phase gas (e.g., air) conditions, with pressure \(P_0\), and with a time-dependent sink/source term given by Eq. (44) that would be specified through tabular generation data in data block GENER. In order for the pressure conditions in this block to be negligibly affected by heat and mass exchange with the flow domain, the volume \(V\) should be made very large, e.g., \(V\) = 1050 m3, just as for time-independent Dirichlet boundary conditions. The time-dependent generation rates from Eq. (44) will then also be very large. The same approach as just outlined for a time-dependent gas pressure boundary can be used to realize other time-dependent Dirichlet conditions, such as prescribed temporal variation of temperature, capillary pressure, and others.
The second, newly implemented method to specify time-dependent Dirichlet boundary conditions is to read a set of time and primary variable data at boundary elements from the input file or a file whose name is given in data block TIMBC. Boundary values will be linearly interpolated between table entries. The volume for these elements should be made very large as well.
Atmospheric Boundary Conditions#
Often flow models need to include the atmosphere as the top boundary. The atmospheric boundary conditions can be specified using Dirichlet boundary elements with very large volumes. For EOS9, no atmospheric boundary element is needed since EOS9 uses Richards’ equation to describe variably saturated flow of a single aqueous phase and treats the gas phase as a passive bystander at constant pressure. A single atmospheric element can be connected to all elements at the ground surface, and users may use AddBound.exe for this purpose (free software available from the TOUGH website).
Atmospheric pressure and temperature are used as initial condition in atmospheric elements. A relative humidity of 100% is conveniently specified by initializing the atmospheric block as a two-phase point with a liquid saturation smaller than the residual liquid saturation (so relative permeability is zero, preventing liquid flow into soil). For relative humidities less than 100%, a single-phase gas point must be specified with an appropriate air mass fraction (1.0 for dry air; the minimum value depends on vapor pressure, which is a function of temperature; intermediate values determine relative humidity). The material properties, such as relative permeability and capillary pressure functions, for the atmospheric boundary element should be appropriately selected so that (1) liquid relative permeability is zero, (2) gas relative permeability is one, and (3) capillary pressure is zero. In addition, mobilities need to be upstream weighted.
Infiltration can be simulated by specifying the infiltration rate in a row of elements below the atmospheric boundary element using the GENER block. Evaporation can be simulated (a) as a binary diffusion process (atmosphere at < 100% relative humidity), (b) by specifying the evaporation rate in a row of elements below the atmospheric boundary element using the GENER block, or (c) by assigning a capillary pressure according to Kelvin’s equation in the atmospheric element (Ghezzehei et al., 2004).
Constant Temperature Boundary Conditions#
Constant temperature boundary conditions but variable pressure, saturation, and mass fractions can be specified by setting rock grain density (DROK) to a very large value (e.g., 1.0E50).
This option is useful for cases where the temperature of the injected fluid is maintained constant while the other primary variables defining the system state may vary with time.