PB91-191692

                                                        EPA/600/2-91/020
                                                                May 1991
                           MOFAT:
      A TWO-DIMENSIONAL FINITE ELEMENT PROGRAM
FOR MULTIPHASE FLOW AND MULTICOMPONENT TRANSPORT

             Program Documentation and User's Guide
                              by

           A. K. Katyal, J. J. Kaluarachchi and J. C. Parker
      Center for Environmental and Hazardous Materials Studies
          Virginia Polytechnic Institute and State University
                 Blacksburg, Virginia 24061-0404

                     Project No. CR814320
                        Project Officers:

                           J. S. Cho
             Processes and Systems Research Division
         Robert S. Kerr Environmental Research Laboratory
                     Ada, Oklahoma 74820
                          L. G. Swaby
               Office of Research and Development
              U. S. Environmental Protection Agency
                     Washington, D.C. 20460
     ROBERT S. KERR ENVIRONMENTAL RESEARCH LABORATORY
             OFFICE OF RESEARCH AND DEVELOPMENT
            U. S. ENVIRONMENTAL PROTECTION AGENCY
                     ADA, OKLAHOMA 74820

                         REPRODUCED BY
                   U.S. DEPARTMENT OF COMMERCE
                       NATIONAL TECHNICAL
                       INFORMATION SERVICE
                       SPRINGFIELD, VA 22161

-------
                                       ABSTRACT
This report describes a two-dimensional finite element model for coupled multiphase flow and
multicomponent transport in planar or radially symmetric vertical sections. Flow and transport of
three fluid phases — water, nonaqueous phase liquid (NAPL) and gas - is considered by the
program which also handles cases in which gas and/or NAPL phases are absent in part or all of
the domain at any given time. The program will simulate flow only or coupled flow and transport.
The flow module can be used to analyze two-phase flow of water and NAPL in a system with gas
present but at constant pressure, or explicit three-phase flow of water, NAPL and gas at variable
pressure. The transport module can handle up to five components which partition among water,
NAPL, gas and solid phases assuming either local equilibrium interphase mass transfer or first-
order kinetically controlled mass transfer. The governing equations are solved using an efficient
upstream-weighted finite  element scheme. Required input for flow analyses consists of initial
conditions, soil hydraulic properties, fluid properties, time integration parameters, boundary
condition data and mesh geometry. Three phase permeability-saturation-capillary pressure
relations are defined by an extension of the van Genuchten model which considers effects of oil
entrapment during periods of water imbibition. For transport analyses, additional input data are
porous media dispersivities, initial water phase concentrations, equilibrium partition coefficients,
component densities, diffusion coefficients,  first-order decay coefficients, mass transfer
coefficients (for nonequilibrium analyses) and boundary condition data. Time-dependent boundary
conditions for the flow analysis may involve user-specified phase heads at nodes or phase fluxes
along a boundary  segment with zero-flux as the default condition. For transport analyses, initial
conditions are specified in terms of equilibrium water phase concentrations of each partitionable
component. Time-dependent boundary conditions may be stipulated as equilibrium water phase
concentrations in the porous medium as prescribed fluxes defined in terms of a specified
concentration in the influent liquid, or with zero dispersive flux specified. Program output consists
of basic information on input parameters,  mesh details and initial conditions plus pressure heads,
saturations and velocities for each phase at every node for specified output intervals. For transport
analyses, the phase concentrations at each node are output at each print-out interval.
                                            11

-------
                                     CONTENTS
       ABSTRACT	ii
       CONTENTS	iii
       FIGURES	v
1.      INTRODUCTION	 1
2.      ANALYSIS OF MULTIPHASE FLOW	2
       2.1 Mathematical model for multiphase flow  	2
        211 Governing equations 	2
        2.1.2 Initial and boundary conditions	3
        2.1.3 Constitutive relations for flow model  	4
       2.2 Numerical model for flow equations	 8
        2.2.1 Finite element formulation	 8
        2.2.2 Treatment of nonlinearity	 11
        2.2.3 Convergence criteria and time-step control	 13
        2.2.4 Adaptive solution domain algorithm	 14
        2.3.5 Elimination of oil phase equation 	 14
3.      ANALYSIS OF MULTICOMPONENT TRANSPORT	 15
       3.1 Mathematical model for transport	 15
        3.1.1 Governing phase transport equations	 15
        3.1.2 Phase-summed equation for local equilibrium transport	 18
        3.1.3 Consideration of nonequilibrium interphase mass transfer	20
        3.1.4 Initial and boundary conditions	22
       3.2 Numerical model description	 24
        3.2.1 General solution approach	 24
        3.2.2 Finite element formulation	 26
        3.2.3 Interphase mass transfer and density updating	 29
4.      GUIDE TO USE OF MOFAT	 31
       4.1 Estimation of soil and fluid properties	 31
        4.1.1 Soil properties 	 31
        4.1.2 Bulk fluid properties 	 32
        4.1.3 Component properties  	 35
       4.2 Program structure and operation 	 38
        4.2.1 Basic features and program control	 38
        4.2.2 Spatial discretization and mesh generation 	 42
        4.2.3 Specification of initial conditions 	 43
        4.2.4 Specification of boundary conditions  	 44
5.      PROGRAM APPLICATIONS	 48
       5.1 Example 1.1-D spill of two-component hydrocarbon	 48
       5.2 Example 2. 2-D planar spill with sloping water table	 50

                                         iii

-------
     5.3 Example 3. 2-D radial dense solvent spill with vacuum extraction	 51
6.                 	 55
APPENDIX A. DESCRIPTION OF MAIN INPUT FILE	 56
APPENDIX B. MAIN DATA FILE FOR EXAMPLE PROBLEM	68
APPENDIX C. MAIN OUTPUT FOR EXAMPLE PROBLEM 	71
APPENDIX D. FLOW CHAJRT FOR PROGRAM	 98
APPENDIX E. DESCRIPTION OF PROGRAM SUBROUTINES	 100
APPENDIX F. DESCRIPTION OF IMPORTANT PROGRAM VAJRIABLES ... 102
                                  IV

-------
                                       FIGURES
Figure 4.1.  Local side numbering of quadrilateral elements	42

Figure 4.2.  Finite element mesh with numbers in parentheses indicating
           global element numbers and others referring to global node  	43
           numbers.

Figure 5.1  Water and total liquid saturation vs depth at end of Stages
           1, 2 and 3 for Example 1	49

Figure 5.2  Aqueous xylene and toluene concentrations at outflow boundary
           for Example 1 for equilibrium and nonequilibrium analyses. Rate
           coefficients in d"1	50

Figure 5.3  Oil saturation contours at the end of Stage 2 for Example 2	51

Figure 5.4  Water phase concentration contours (gm~3) at the end of
           Stage 2 for Example 2	51

Figure 5.5  Gas phase concentration contours (g m'3) at the end
           of Stage 2 for Example 2	51

Figure 5.6  Oil saturation contours at the end of Stage 2 for Example 3	53

Figure 5.7  Water saturation contours at the end of Stage 3 for Example 3	  53

Figure 5.8  PER mass in system vs time during Stage 3 of Example 3	53
                                           v

-------
                                   1. INTRODUCTION
This manual describes a two-dimensional finite element model, MOFAT, for coupled multiphase
flow and multicomponent transport in planar or radially symmetric vertical sections. Flow and
transport of three fluid phases—water, nonaqueous phase liquid (NAPL) and gas—is considered
by the program which also handles cases in which gas  and/or NAPL phases are absent in part or
all of the domain at any given time.

The program allows the user to analyze flow only or coupled flow and transport. The flow module
can be used to analyze two phase flow of water and NAPL in a system with gas present but at
constant pressure, or explicit three phase flow of water, NAPL  and gas at variable pressure. The
transport module can handle up to five components which partition among water, NAPL, gas and
solid phases assuming either local equilibrium interphase mass transfer or first-order kinetically
controlled mass transfer. Time-lagged interphase mass transfer rates and phase densities are used
to avoid costly iterative procedures between the flow and transport analyses while enabling
accurate mass balance and providing for the effects of density-driven flow.

The systems of governing equations for multiphase flow and multicomponent transport are solved
using an efficient upstream-weighted finite element scheme. To avoid costly numerical
integration, an influence coefficient approach is utilized using linear rectangular elements.
Nonlinear time integration in the flow analysis is handled using a Newton-Raphson method with
an implicit saturation derivative formulation of the governing equations which is very efficient and
accurate.

Required input for flow analyses consist of initial conditions, soil hydraulic properties,  fluid
properties, time integration parameters, boundary condition data and mesh geometry. Three phase
permeability-saturation-capillary pressure relations are  defined by an extension of the
van Genuchten model which considers effects of oil entrapment during periods of water
imbibition. The model requires specification of parameters defining the air-water capillary
retention function, NAPL surface tension and interfacial tension with water, NAPL viscosity,
NAPL tensity, maximum residual NAPL saturation and soil hydraulic conductivity.  The latter may
be anisotropic  and soil properties may vary spatially.

For transport analyses, additional input data are porous  media dispersivities, initial water phase
concentrations, equilibrium partition coefficients, component densities, diffusion coefficients,
first-order decay coefficients, mass transfer coefficients (for nonequilibrium analyses) and
boundary condition data.

Initial conditions for the flow analysis are described by phase pressures expressed in terms of
water-equivalent heights. If the analysis involves NAPL, the program has the capability to
determine and assign NAPL heads which correspond to the zero oil saturation condition. Time-
dependent boundary conditions for the flow analysis may involve user-specified phase heads at
nodes (type-1 condition) or phase fluxes along a boundary segment (type-2 condition) with zero-

-------
flux as the default condition. For transport analyses, initial conditions are specified in terms of
equilibrium water phase concentrations of each partitienable component. Time-dependent
boundary conditions may be stipulated as equilibrium water phase concentrations in the porous
medium (type-1 condition) or as prescribed fluxes defined in terms of a specified concentration in
the influent liquid (type-3 condition). The program assumes zero dispersive flux (type-2
condition) along any boundary segment not explicitly specified.

Program output consists of basic information on input parameters, mesh details and initial
conditions plus pressure heads, saturations and velocities for each  phase at  every node for
specified output intervals. In addition, the total volume or mass of each phase, time-step size and
number of iterations are printed at each time-step. For transport analyses, the phase concentrations
at each node are output at each print-out interval.  The program allows the user to restart
simulations from a previous run using an auxiliary data file automatically dumped from the
previous simulation to define initial conditions for the restart problem.

Chapters 2 and 3 of this document describe the model development and numerical  methods
pertaining to multiphase flow and multicomponent transport, respectively. Chapter 4 discusses the
estimation of soil and fluid properties required for the model and describes various practical
aspects of program usage. Detailed instructions for creating the data file needed to run the
program, and example input and output files are given in appendices.

-------
                         2. ANALYSIS OF MULTIPHASE FLOW

2.1 Mathematical Model for Multiphase Flow

2.1.1 Governing equations

The mass conservation equations for water (w), organic liquid (o) and air (a), assuming an
incompressible porous medium, incompressible liquid phases and compressible gas phase, may be
written in summation convention for a two dimensional Cartesian domain as (Parker, 1989)
               *
                 as.     9f..   *.                                                  [21a]
                  dt      dx±
                                   Ro
                  dt      dx1     Po
                   dt         dx^      a



where (() is porosity, 5^ is the p- phase saturation, x;(and Xj) are Cartesian spatial coordinates (i,j =
1,2), g  is the Darcy velocity of phase p in the/'-direction, ppis the density of phase p, R^is the
net mass"transfer per unit porous media volume into (+) or out of (-) phase p, and t is the time.

Darcy velocities in the p-phase are defined by


                            J d h           }
                 "p. ~    P. . 1 3 .,  + "rp  if                                         \~-^\
where Kp   is the p-phase conductivity tensor, hp = Pp /gp"w is the water height-equivalent
pressure nead of phase p where Pp is the p-phase pressure, g is gravitational  acceleration and p
is the density of pure water, pp is the density of phase p, pr  =  p  /p*K is the p-phase specific
gravity, and u^ = dz /dxj  is a unit gravitational vector measured positive upwards wherez is
elevation.

Combining [2.1] and [2.2] yields

-------
                  at    ax
          K
           "w« lax.
                                                                                   [2.3a]
ss0 _   d
 at    ax.
                                      ah
                                                                                   [2.3b]
                 at
     ax.
                                                                                   [2.3c]
If fluid saturations are treated as implicit functions of the phase pressures, equations [2.3] may be
solved directly with phase pressures as primary variables. Details of the numerical methodology
will be described in a subsequent section.

The system of flow equations solved may be reduced at the user's discretion. In many cases, gas
pressure gradients may be regarded as negligible and the gas flow equation may be discarded. In
this case, a gas phase may still be present and may be mobile, but it is assumed that gas flow has
negligible effect on liquid flow due to small gas pressure gradients (Richard's assumption). The
flow equations may be written also for the case of a two-dimensional radially symmetric domain
using standard coordinate transformations. MOTRANS provides the capability to perform
calculations for either Cartesian or radially symmetric domains.  Options in the flow module can
thus accommodate analyses for vertical plane slices (Cartesian symmetry) or vertical radial slices
with two options for stipulation of the flow equations as follows:

1) Three phase flow of water, oil and gas phases (gas pressure variable)
2) Two phase liquid flow of water and oil (gas phase at atmospheric pressure)

Details of model  development will be presented in this report for Case 1 with Cartesian geometry.
Case 2 has an identical form to that of Case 1, but is degenerate in that the gas flow equation is
omitted. The analysis of radial flow is likewise similar except for minor changes in the governing
equations in converting them from Cartesian to radial symmetry.

2.1.2 Initial and boundary conditions

Initial conditions for each phase p must be stipulated on the entire flow domain,/? as
            hp(x.,0)  - hpl(x.)
                     on  R  for   t = 0
                                                                                    [2.4]

-------
Boundary conditions may be stipulated as type-1 or type-2 as follows
          h  ( x , t)  = h  (x  , t)       on  S   for  t>0
           p    i        p2  *                 l                                    [2.5a]

             *Plni = ^Pi(xi' tj        on  S2  for   t>0                          [2.5b]
Equation [2.5a] denotes type-1 boundary conditions (pressure head) on boundary region $. The
type-1 condition will commonly be used for the water phase on boundaries within a saturated
aquifer where positive heads are specified corresponding to assumed vertical hydrostatic
conditions (i.e.., hw = zwt - z where zw,is the water table elevation and hwis the water head at
elevation z). Type-1 conditions may also be used on infiltrating boundaries for water or oil if the
infiltration is pressure controlled (e.g., surface spill with free liquid on surface). Gas phase
boundary conditions will often be type-1 corresponding to zero pressure along surfaces in contact
with the atmosphere or at negative or positive heads at vacuum or injection points.

Equation [2.5b] denotes type-2 (fluid flux) boundary conditions with qp (x.,t)  representing the
normal flux on boundary segment S2 for phase p with outward normal unit vector nt. The default
flow boundary condition is zero  normal  flux. Nonzero flux boundary conditions will be commonly
used on boundary regions along  which inflow rates are known. Specifying fluid removal rates is
also possible, but may be somewhat troublesome since removal rates are often controlled by soil
permeability to the phase which  is not known a priori.

MO TRANS permits stipulation of type-1 or type-2 boundary conditions for all fluid phases
independently using a very flexible scheme which utilizes user-defined schedules to specify time-
dependent behavior. For problems run piecewise, a restart option is available.  Mixed type
boundary conditions can occur in certain problems, but are not considered directly. Such
circumstances can arise,  for example, when water or oil seeps from a boundary which is open to
the atmosphere. The correct boundary conditions in such a case would involve stipulation of zero
fluid flux as long as the boundary head remains  negative, but switching to a type-1 condition with
zero pressure after the pressure build up to atmospheric. To model such circumstances, the user
will need to stop the code and restart it after manually changing the boundary  conditions.

2.7.3 Constitutive relations for flow model

Relationships between phase permeabilities, saturations and pressures are described by a three
phase extension of the van Genuchten model which takes into account effects of NAPL entrapment
on NAPL relative permeability. To describe the saturation-capillary pressure relations, we
introduce capillary pressure heads defined by

                              h« - h* ~ h«                                          [2.6a]

                              h°"  = h°~h«                                          [2.6b]

-------
                                                                                      [2.6c]
where hp = P/p^g with Pp the/>-phase pressure, pw the density of water, and g the gravitational
acceleration andp = a (air), o (oil) orw (water).

Prior to the occurrence of oil at a given location, the system is treated as a two-phase air-water
system described by the van Genuchten (1980) function
where sw =  ( sw -  sffl) / ( 1  - sm)  is the "effective" water saturation, S^is the "irreducible" water
saturation, a [L"1] and n[-] are porous medium parameters and m=\-\/n.

Following the occurrence of oil at a location, the system is described by the three-phase relations


                       S   = [ 1 +  ( a P   h   ) n rm
                        W   L     V   ~ O W  nitf '  J
                                      ow  ow'
st =
                                                                                      [2.8a]

                                                                                      [2.8b]
where s w is the apparent water saturation, s t is the effective total liquid saturation, Pao is a
scaling coefficient which may be approximated by the ratio of water surface tension to oil surface
tension and  P^ is a scaling factor approximated by the ratio of water surface tension to oil-water
interfacial tension. The apparent water saturation is given by
                              sw = Sw + sot                                            [2.9]
where sot is the effective trapped oil saturation. The latter is estimated using an empirical
relationship given by Land (1968) as
                i- sm.in              i- s.
    Sot- K  _   ..  ^in.     .       ..  =  . I    	V  -.                    [210a]
                                                                                     [2.10b]

-------
where

                            Ro«= -^r ~  l                                        p.ioc]
where s™*x  is the maximum effective residual oil saturation where effective oil saturations are
expressed in the form S0/(l-SJ. The minimum effective water saturation, ~s*in , is defined as the
historical minimum effective water saturation at a given location since changing from a two phase
air-water to a three phase air-oil-water system. The value of ~s*in  is updated at the end of each
time step after convergence and stored for each node. Actual water and oil saturations may be
computed as
                                        - sot)  + sm
and

                      So= (1  - Sm>  St+ SK-  S»                                  [2_nb]
To avoid sudden discontinuities at a node due to the change from two phase air-water to three
phase air-water-oil system under condition for which l/pao+l/pow *1, a node once classified as a
three phase node will not be allowed to be considered as a two phase air-water node (Kaluarachchi
and Parker,  1989). The necessary condition to classify a given node as a three phase node is
                                                                                   [2.12a]
where
                                        3
h
 o
                           cr    30  a   ^ ow  w

                                   POW+Pao                                        [2.12b]
Phase conductivities are described by
                            K    = k_  K
                                                                                   [2.13a]

-------
                         x\    — K.   x\     / ll                                        r^ I o 1 ~\
                           o±.    ro  sw..  /  Iro                                     [2.13b]
                         K   = k   K     / T\
                           a      r a  aw. .  /  '
                           a±j    ra
«i, /  '«                                      [2.13c]
where krp is the relative permeability of phase p, TJ^ is the absolute viscosity ratio between phase p
and water, and K     is the saturated conductivity tensor for water. It is assumed that the
                 j. j
coordinate system is oriented with the conductivity tensor, or otherwise that off-diagonal
components may be disregarded, so that KSW  = 0 for i * j.

Relative permeability relationships corresponding to the three phase van Genuchten model have
been derived by Parker et al. (1987).  These may be modified to account for effects of NAPL
entrapment as described by Kaluarachchi and Parker (1991). The relative permeability relations
are given by

                     lr   - ~C ^ I" 1   M   ~Q 1^m \ m 1 2
                       rw ~wL'w'J                                   FO1/I1
                       rw     w             w                                         [2.14aJ

                 —    =   %        = 1/m  m       — 1/ffl  m 2
         kro=(st-s«>    Ld-sv    )   -  (i-st   )   J                           [2.14b]
                                          •^ 1/m , 2m
                                            c                                       [2.14c]
where parameters are all as previously defined for the saturation-pressure relations.

When gas phase flow is considered, the gas compressibility is described by the linear relation

                           p  = Ah  + p° + pc
                                   a  ^a   ^a                                        [215]

where A is the gas compressibility taken to be 1.17 x 10"6 g cm'4, pa° is the density of native soil air
taken to be 1.12 x 10"3 g cm'3, and pac is the density of contaminants in the gas phase which is
determined from the solution to the transport model when coupled flow and transport is analyzed
(see Chapter 3).

-------
2.2 Numerical Model for Flow Equations

In this section, we describe the numerical model for the flow equations. The model description will
be given for the case of fully three phase flow in a two-dimensional plane vertical slice. Analogous
procedures are applied for the case of radial flow domains and for two phase flow with gas at
atmospheric pressure.

2.2.7 Finite element formulation

The governing equations in this formulation have been solved using a finite element method with
spatial derivative terms in the flow equations approximated by asymmetric upstream-weighting
functions developed by Huyakorn andNilkuha (1978), while the remaining terms are handled
using linear basis functions. For efficiency, the model is restricted to the use of linear quadrilateral
elements. The solutions for pressure heads over an element are approximated by trial functions of
the form
                               Y,  wTU,ii) hwj(t)                              [2.16a]
                               1=1
                  ho (x., t) - £,  NT(^,T\) hol(t)                              [2.i6b]
                  ha(*lft) - 1,  tfTU,H) haj(t)                              [2.16c]


where hwl, hol and hal are pressure heads for water, oil and air, respectively, at node/, N, is the
linear basis function for node /, and £ and TJ are the local coordinates. Applying Green's theorem
and the standard Galerkin principle to [2.3] with upstream-weighting for spatial terms yields

                      dW, dN,          ,    dN,
             N, Nj      dR +    qnw N, dS - 0    /, J- 1, 2,... N                     P-17a]
                     dW   dN                 dN
                 K		 h T dR + IK	 u . dR
               dS
               —^ dR +/  
-------
                              dW   dN                     dN
                    I Pa Ka  T— ^  T— ^ haJ dR +  I  Pa Ka  ^  Ui d*               [2.17c]
                    J  a   ai  dx.  dx.   aj       J   a  ai  ax.   1                  L     J
                              dx.  dx.
                    R
                     d(p  s
                                 dR + J gna Wr dS = 0     J,  J= 1,  2, .  .  . N

where fFj is an upstream-weighted shape function described in detail by Kaluarachchi and Parker
(1989) and N is the number of nodes in the domain. Since only diagonal terms in the conductivity
tensor are considered in the present analysis, the double subscript notation is henceforth dropped.

Following Kaluarachchi and Parker (1989), we expand the matrices in the first two integrals of
[2.17a], [2.17b] and [2.17c] into four submatrices by permitting the hydraulic conductivity within
an element to vary through the upstream-weighting function such that
                   *„  («,., t)  = x,  wx u,n)  V
                                                                                    [2.18]
                    p<    *       r.i
where K  denotes the hydraulic conductivity attributed to node/. In the present analysis, we take
K r v
gives
K r values corresponding exactly to node points. Substituting [2.18] into [2.17] and integrating
                                 K
              A T T h  T + BT T 	— - F T =  0
                w IJ  w J    IJ    I-     w I
                                                                                   [2.19a]
                                    dS
                                    -7^-FoJ-0                                 [2.19b]
                                     d t
                                     dt

-------
where
                    /
                         4                ,4

                                    W"   m
                     ^ E   [tf^j/^-
                                  J  '      g.1           y                      [22Qa]


                                  E  V                                      P.a«]
                      n    ,  4
                wj = - E  - E   lGleT K1 + boundary flux term                [2.20e]
                     e=i   2 r=i          y
                      n    ,      4
               F  = - E  - Pro E  [G] r Ko + Boundary flux term             [2.20f]
                     e=l   2     r=i          y
                      n   d      4
                  = - E  — PrE  [G]jK'I+ boundary flux term            [2.20g]
                     e = l   2   r« 1 = 1         ay
where w and J are the width and length of the rectangular element respectively. Influence
coefficient matrices,  [ u*]^,  [ uy]eT,  [ G ] ^ and coefficient vectors for boundary flux terms are
given by Kaluarachchi and Parker (1989).  In the subsequent analysis the boundary flux terms in
[2.20e] - [2.20g] have been omitted for simplicity. A mass lumping procedure is employed to
diagonalize the mass storage matrix.

                                          10

-------
2.2.2 Treatment of nonlinearity





A Newton-Raphson iteration method is used with this formulation to solve the nonlinear problem.


In residual form, [2.19] may be written as



                                         dS j

                        ~-        «j + Bu —TT-  - FWI

                                                                                    [221a]
                                          dS  7

                                         --FOI                                P.21b]
                                           dt
Applying the Newton-Raphson procedure to [2.21] gives
                          >1    r     9rwl     r.l
                 ™.j    ^     ^    9hoJ
                                                       T

                                                     °J
                   oJ/,r*l   i,rv       or
                        dh  T
                          a J
                                    -
                                 a J    a J
                          r+1   i,-^\       a-^/i,-f+l   i,-f\
                           -r  -  h  T)  +  	  (h  T  - h  T)
                          'J    w J'    a u    *oJ     oJ'
                   wJ                    oJ




                           3 -i  /  j_ JT+1   i_ JT \
                        dh
                                                                                    P.21c]
                                         OJ



                          :DI "-^1   - r ,    _ rr                                  [2.22b]
where r and r+1 refer to the previous and current iterations, respectively. The derivatives in [2.22]


can be expressed as
                                             11

-------
                h  . 	— + [BTJ 	— ] /At	—
                 W Ll  *J L.        -L (J 3 T_           3 Ji

                     5h^         5h^         3h"*                 [2.23a]
                             as  T
                            -^?] /A
               dA TT         as  T
         = h    	°J± + [B   	££] /At
                                12
                                                                     [2 23b]
                                   /At-—^                       [2.23c]
                                                                     [2.23d]
                                as
       _      ,      O ± Li   I-TI      O <-* n  / A j_
     = A    + h    	  + [ B    	 ]  / A t -
3L.      oJJ    o-^3j-,        IJ 3 i^            3Z-,
on  T               on   T         0/2  T          on T                   .... _-  ,
  oJ                 oJ           oJ            oJ                   ^2 23G!
                                                                     [223g]
                 arr           ^i            aT
   	i£  = h  r 	iit + [BrT 	^ ] /At	^
   ah T    aL  ah  T     JJ ah  T         ah T
      °J          °J           °J            °J                       [2 23 hi

-------
dr  T
  3 I    ,
      = A
    T
  aj
                                    dA  Tr          dS  T          dF T
                               T_      a I L    r _,      <3 J .,  / A ,       B I
                             + h   -  +  [ B   - ]  / A t - -
                                aL  dh   T      IJ  dh  7          dh T
                                       aj            aj             aj               [2 2311
where At = f+l - f and the subscript L is a dummy index denoting the appropriate influence
coefficient matrix to be selected for a particular set of /and J. In [2.21],  dswj,  dsoj and  dsaj
are changes in the water, oil and air saturations, respectively, over the time interval calculated as


                            d S   = S *~*^ -  S *"
                              •TJ-  wj     wj                                        [2.24a]
                                 ~~
                                                                                       24b]
                            dSaJ= SaJ  ~  SaJ                                        [2.24c]

where t and ^+ 1 denote previous and current time-steps. Phase saturations are calculated from
phase pressures using the constitutive relations described earlier. Saturations at the t+ 1 time-step
in [2.24] as well as residuals rwl and rol in [2.21] and derivatives with respect to pressure heads in
[2.23] are evaluated using currently updated pressures for the present time-step. Equation [2.22] is
then solved for the changes in pressure heads, the nodal pressure heads are updated, and the
procedure is iterated until the solution converges.

2.2.3 Convergence criteria and time-step control

The convergence criterion for phase p (= o, w, a) is as follows

                       I hrtl - h r\ <.  I h r\ e  + e
                       IP     Pi   i  p  i   r    a                                    [2.25]

where r+ 1 and r refer to current and previous iterations, er is the allowable relative convergence
error and ea is the allowable absolute convergence error for phase p. The convergence criteria must
be met for all phases and at all nodes.

Automatic adjustment of the time-step size is undertaken depending on the number of iterations
required for nonlinear convergence relative to a specified maximum number of iterations, Imax,
given by the user. If the number of iterations is less than or equal to!max/4 the time-step size is
increased by a factor F for the next time-step. If the number of iterations exceeds lmjl> but is less
than Imax, the solution proceeds to the next time step but the time step size is decreased by a factor
F. If the number of iterations exceeds Imax, program execution is halted and output is dumped to
enable the problem to be restarted (e.g., after reducing the time-step size). For most simulations,
good results may be obtained using/max.= 30 and F= 1.03 to 1.05.

                                             13

-------
2.2.4 Adaptive solution domain algorithm

In many practical problems involving soil contamination with immiscible organic liquids, a
substantial part of the spatial domain may exhibit negligible changes in fluid saturations over all or
part of the time domain. For example, during infiltration of oil, the region well ahead of an oil
front may be essentially static. Also, after redistribution proceeds to the point "residual saturations"
are  achieved or under steady flow boundary conditions, quasi-static or steady-state conditions may
occur which lead to subregions with negligible changes in fluid pressures and saturations. Solution
of the flow equations in such subregions involves unnecessary computational effort. An adaptive
solution domain  (ASD) formulation was developed to confine the mathematical solution domain to
a subdomain within which transient oil flow occurs in order to improve computational efficiency.
The method involves reducing the solution domain for the water and oil flow equations by
eliminating from the global matrix assembly elements which are not experiencing significant
changes in fluid  saturations and/or pressures. The formulation achieves computational efficiency
by solving the flow equations only for the part of the domain where changes in fluid pressure or
the  saturations take place above stipulated tolerances.

In the ASD method, elements are categorized as either "active" or "inactive" at a given iteration. If
the  element is classified as active, then it is included in the global matrix  assembly, otherwise it is
excluded. The criteria for changing a node from inactive to active is based on the changes in water
pressure head, oil pressure head, oil saturation, and water saturation at connected nodes from the
last converged time-step to the current iteration. The  specific criteria for changing an element from
inactive to active is as follows
                    > TOL1                                                      [2.26a]
           or
               A/z0 > TOL1                                                      [2.26b]
           or
               &SW > TOL2                                                      [2.26c]
          or
               &S0>TOL2                                                       [2.26d]
whereA/7w, A/70, A5"w and A 5^ are the changes in hw, h0, Sw and S^ respectively, between the
previous converged time-step and the current iteration, TOL1 is a specified head tolerance for the
ASD method and TOL2 is a specified saturation tolerance. Equations for nodes on the active-
inactive frontier will have contributions only from active elements and the boundary will shift
naturally as imbibition or drainage "fronts" pass through the domain.

2.3.5 Elimination of oil phase equation

In many practical situations involving flow of organic liquids, the oil phase fluxes after long
redistribution periods become negligible, resulting in insignificant changes in oil phase saturations

                                            14

-------
within the entire domain. MOTRANS has an option to globally eliminate the oil equation from the
solution if changes in oil saturation at all nodes during the current time step are smaller than
tolerance TOL3.
                                            15

-------
                  3. ANALYSIS OF MULTICOMPONENT TRANSPORT


3.1 Mathematical Model for Transport

3.7.7 Governing phase transport equations

In this chapter, we will describe a model for transport of NAPL constituents which can partition
among water, oil, air and solid phases following/'arfer (1989) and Kaluarachchi and Parker (1990).
To model component transport, continuity and mass flux equations for each partitionable component
in each phase must be specified. Mass conservation of species a in the p-phase requires that
                        a ..         a        ap    ap                               n 11
                        3 t         3 x^                                            p-lj


where Cap is the concentration of the noninert a component in p-phase expressed as the mass of a
per phase volume [M L"3], Jttpi is the mass flux density of a in p-phase per porous media cross
section in the i-direction [M L"2 T1], Rap is the net mass transfer rate per porous medium volume of
species a into (+) or out of (-) the p-phase [M L"3], and yap is the net production (+) or decay (-) of
a within phase p per porous medium volume due to reactions within the p-phase [M L"3] described
subsequently by


                             ' ap ~ ~ "ap ap                                         |-~ ~-,

where nop is an apparent first-order decay coefficient.

The mass flux density of component a in phase p due to convection, diffusion and mechanical
dispersion is described by


                                                aCap
                                                gxj                                [3^



where Dapij is a dispersion tensor given by

                      n       ~   I  r^dif   r^hvd \
                      D   . .  = Q   (  D     + LJ .  . 1
                       dpi]   -'ap   ap     pi]                                    f'J A~\

in which Ddif is the molecular diffusion coefficient of a in the p-phase of the porous medium, D

                                           15

-------
is a mechanical dispersion coefficient, andg^ is a nondilute solution correction factor.  For the case
of transport of low solubility organic components, small volume fractions of organic components
in water and gaseous phases will occur and gaw =  gOLa = 1 is assumed.  For nonaqueous phase
liquids, the phase composition may reach 100% (e.g., single component organic liquid) at which
point dispersive transport becomes nonexistent.  We accommodate nondilute solution diffusion-
dispersion by approximating the oil phase nonideal solution factor by
                                        P°                                          [3.5]
Employing the tortuosity model ofMillington and Quirk (1959) yields the expression for o£lf as
                          Ddlf - d>1/3  S7/3  D°
                           ap    ^     p    ap
where D° is the diffusion coefficient of a in bulk p-phase.  The mechanical dispersion
                                                                                    [3.6]

                                                 )hase.  The mechanical dispersion
coefficient is shown by Bear (1972) to have the form



                                                 LT>    qp-qpj                       [3.7]
                       T -" n                              Q
                          tr                               ]P



where AL and^4r are longitudinal and transverse dispersivities [L], qpi and qpi are p-phase Darcy
velocities in the i and j directions,  g  =  Eg2.]1* is the absolute magnitude of the p-phase
                                 p
velocity, and 6tj is Kroonecker's delta.

Combining the phase continuity equation [3.1] and the mass flux equation [3.3] for transport of
component a in the p-phase yields


    .    CapSp    d                 Cap       Cqp
         dt      dx       P
Expanding the first and third terms in [3.8], employing the bulk p-phase continuity equation
[2.1], and assuming density derivative terms to be of second order importance within a given
time step yields
                                           16

-------
4) sr
      ac
        ap
       at
   ax,
              S D
               p o
                                    ax,
                                                                     ap
                                                                          [3.9]
  where the total phase mass transfer rate, Rp, it may be noted, is related to the individual
  component mass transfer rates by

                               *P •  £  v
                                                                                   [3.10]
  where n, denotes the number of "noninert" or partitionable species. In the present context, we
  will use the term "inert" to refer to components of the NAPL phase which are for practical
  purposes insoluble and nonvolatile. For the three fluid phase system, [3.9] constitutes a system
  of three equations which we may write for the water phase (p=w), the organic liquid phase (p=o)
  and the gas phase (p=a) as
         dc
"  at     ax.
                                 ax.
                                          4wi
                                    ac
                                      aw
                                    ax.
                                                           R
   «D
at     ax.
                                              dx.
                                                            H1 fy f.
                                                                                  [3.lib]
   at     ax.
                      C(>S aD
acM
 ax.
                                              ax
                                                                        [3.lie]
  To accommodate adsorption of a by the solid phase, an additional continuity equation is required
  which may be written as
                            a t
                                    as  "as  as
                                                                         [3.12]
  where Cas is the solid phase concentration expressed as mass of adsorbed component a per
  porous medium volume [M L"3],|ias, is the first-order decay coefficient (T1) in the solid phase and
  Ras, is the mass transfer rate per porous medium volume to (+) or from (-) the solid phase [M L"3
  T4].
                                            17

-------
3.1.2 Phase-summed equations for local equilibrium transport

Coupling between the phase transport equations arises due to the interphase transfer terms.
Explicit consideration of interphase transfer kinetics may often be justifiably avoided by
assuming phase transfer to be equilibrium controlled. We consider first the case of linear
partitioning and introduce the approximate thermodynamic relations
              C
               ce
              C
              ^cs
              c
r   c
 ceo   ccw
r   c
L csa ^ aw
r   c
                                        [3.13a]
                                        [3.13b]
                                        [3.13c]
where rKO is the equilibrium partition cofficient for species a between water and organic liquid
(Raoult's constant),  raa is the equilibrium partition coefficient between water and gas (Henry's
constant), and rKS is a dimensionless equilibrium partition coefficient between water and solid
phase. Using the equilibrium relations, we may rewrite the phase transport equations in terms of
a single phase concentration.  For water-wet systems, it is logical to retain the water phase
concentration, since water will always be present in the system. Using [3.13] to eliminate oil, gas
and solid phase concentrations from [3.11] and summing the equations noting that

             R+R+R+R=0                                            [314]
                tKW    (XO    (%S    CHQ                                                L    J

leads to the a-componentphase-summed transport equation in terms of water phase concentrations
                    d t
                                  :>*..
                                   aij
ax.
                                                                                   [3.15a]
where
                aij   ™  w  awij  ™  o  awij  OLO + ™   a  awij  aa
                                           [3.15b]


                                           [3.15c]


                                           [3.15d]
       *   ,1
         =U
       o<   K
                                                                                  [3.15e]
                                            18

-------
Note that [3.15] has the same form as the simple single phase transport equation. However, the
coefficients represent pooled effects of transport in all phases.  Note also that interphase mass
transfer terms occur in the phase summed equation only as a sum over all components.
An alternative form of the phase-summed transport equation may be written in terms of the oil
phase species concentrations. The a-componentphase-summed transport equation in terms of
oil phase concentrations has the form
               01    at     ax.
                              .
                                D**.
 ** " '"ao
•ai   ax,
                                                                                 [3.16a]


where
                                                  ras
                                r        r        r
                                 ao       ao       ao                            [3.16b]
                   ™
              aij  -1-  o  aoij      r              r                              n  16H
                                    ao            aa                            [J. 1OCJ
                        ## =     +  ^"^    ^ai
                       "oii ~ ^oi   r1
                                           r
                                                                                 [3.16d]
                        i]     i]   T      i]  I"1      /?       R
             ## _       ^aw   ^aa aa    ^as as     o       w
            rTV    rTV C
                        r       r          r       p     p r
                         ao      aa         aa     'o    *» ao
                                                                                 [3.16e]
The form of the oil-based equation is identical to that of the water-based phase-summed
equation. The water-based equation has the conceptual advantage that Cm is physically
meaningful at all locations in the porous medium, whereas Cao only has physical meaning at
locations where S0>0.  The latter constraint, however, imposes no mathematical difficulty since
[3.13] can in any event be employed as a definition of Cao as well as Caa whether or not the phase
exists.  MOTRANS includes provisions for solution of either the water- or oil-based phase-
summed transport equations. Since most of the species mass commonly resides in the oil phase,
the oil-based equation has the advantage that mass balance accuracy of the solution may be more
readily controlled.

                                          19

-------
When simulating two phase liquid flow consisting of water and oil with constant gas pressure,
only diffusive transport is considered in the gas phase— i.e., qa = 0 is assumed for this case.  It
should be noted that even if gas pressure gradients have negligible effect on liquid phase flow,
this does not strictly imply zero gas flow. If fluid saturations change in time, gas  flow will occur.
Furthermore, phase partitioning may result in density gradients in the gas phase which can induce
significant flow.
3.1.3 Consideration ofnonequilibrium interphase mass transfer

The phase-summed formulation of the transport equations may be generalized to account for
nonequilibrium phase partitioning by introducing the concept of apparent partition coefficients.
Under transient field conditions, at a given location in time and space, actual phase concentration
ratios may differ from the true equilibrium ratios defined by [3.13]. With this in mind, we define
apparent partition coefficients — which will vary in time and space — by analogs of [3.13] as
                             Cao = rao Ca«                                        [3.17a]
                                                                                  [3.17b]
                             c  = r'  c
                               as     as  «„                                        [3.17c]

For any two phases which are in physical contact, the rate of mass transfer will be described by first
order mass transfer functions of the form
                        p    -  V    i r     r1 e \
                         0(12 =   «12  (  al "  al '                                   [3.18a]
                                                                                  [3.18b]
where RKl2 is the rate of mass transfer of a per porous medium volume from phase 1 to phase 2,
Cal is the actual concentration of a in phase 1, c^ is the concentration of which would occur in
phase 1 if it were in equilibrium with phase 2, Ca2 is the actual concentration of a in phase 2,
Cae2 is the concentration of a that would occur in phase 2 if it were in equilibrium with phase 1,
and kal2 is a mass transfer rate coefficient [T1].
                                           20

-------
Consider first, the case of mass transfer between oil and water phases when both phases exist at a
point in time and space. Employing [3 . 1 8a] along with [3 . 1 7a] for the actual oil phase
concentration and [3. 13 a] for the equilibrium concentration gives
                                           r   c)
                                            ao   aw'                                 ["3 inn
which may be solved in terms of r'   as
                                         aow
                          ao    ao   7,    r-
                                      aow  aw
                                                                                    [320]
which indicates that the apparent partition coefficient my be expressed as the actual partition
coefficient plus a correction term which depends on the sign and magnitude of the actual mass
transfer rate and on the concentration in the water phase.

If free oil is present in the system at a given point in time and space, so that air and oil are
physically in contact, mass transfer between oil and gas phases may be described in a similar
fashion. Employing in this case [3.18b] we obtain

                   R    = - k     (r'   c   - r   c  )
                    aoa     aoa  v  aa  aw    aa  aw'                              r-j oil
which yields
                                        Raoa
                                                                                    [3 22]
In the absence of free oil in the porous medium at a given location, mass transfer will occur
between water and gas phases rather than between oil and gas phases.  Proceeding in the same
manner as for oil-gas mass transfer yields
                                        Rawa
                                           21

-------
Note that either [3.22] or [3.23] will be relevant at a given location and time, depending on
whether free oil occurs.  Finally, mass transfer between water and solid phase may be considered
by employing [3.18b] to obtain
                          p/  _ p         OfWS
                           as    as                                                 n o/in
                                             aw                                     [3.24]
In order to describe transport with nonequilibrium phase partitioning, the individual phase
transport equations, [3.1 la-c] and [3.12], may be summed as described in section 3.1.2 but using
apparent rather than equilibrium partition coefficients. Thus [3.15] will hold for the
nonequilibrium case with apparent partition coefficients used in lieu  of the equilibrium
coefficients. Nonlinearity is introduced into the resulting phase-summed transport equation due
to the dependence of apparent partition coefficients on concentration as well as mass transfer
rates. Thus, an iterative solution procedure will be required for the nonequilibrium problem as
will be discussed later.

If the transport equations are written in terms of oil phase concentrations; similar expressions to
those given above may be derived to define apparent partition coefficients as functions of mass
transfer rates and oil phase concentrations.
3.1.4 Initial and boundary conditions

Since the transport equation is written in the phase-summed form, it is necessary to specify initial
and boundary conditions in terms of a single concentration.  In the following, as well as in the
program itself, we will specify initial conditions and type-1 boundary conditions in terms of
water phase concentrations. If the oil concentration form of the phase-summed transport
equation is employed, oil phase concentrations are internally computed from the specified water
phase values and the equilibrium relations.  Type-3 or flux-type boundary conditions are defined
in terms of phase-summed fluxes. Initial conditions for component a are given as

                C   (x  ,  t=0 )  =  C    ( x )            on  R  for t  = 0
                 aw    1            a"-    '                                          [3.25a]

where caw  ( x_. )  represents the initial water phase concentration of component a at location xt.
Note that the concentrations of a in all other phases are implicit in [3.25a] via the phase partition
relations [3.13] or [3.17]. Boundary conditions for the phase-summed equation are stipulated by
                KK (xit  t)  - CaK (x±,t)           on S1 for t > 0              [3.25b]
                                            22

-------
                a c
                - —  = 0                     on S2 for  t > 0                     [3.25c]
                                           3 C    /
                     Crire = q* C   -  D*. .  - —   n      on S,  for t >  0
                       ctw     ^ n  ctw    0£_z.j   3    1             3                   ro
Equation [3.25b] describes a type-1 boundary condition on boundary region^ with specified a-
species concentration Cawl (xp t) within the porous medium at the boundary. Note again that
stipulation of the water phase concentration simultaneously fixes the organic liquid, gaseous and
solid phase concentrations via the partition relations. Type-1 boundary conditions for transport
will be used rarely since direct control-over porous media concentration seldom occurs. One
case which may be relevant would be the occurrence of a source of known concentration where
diffusion is the dominant transport mechanism (a type-3 condition would be used if significant
convection occurred at the boundary).

Equation [3.25c] describes type-2 boundaries with zero normal concentration gradient specified
on region S2. The zero gradient condition indicates zero normal dispersive flux at the boundary.
If normal fluid velocities are nonzero, transport across the boundary will be permitted by
convection.  This boundary condition is commonly applied at boundaries which have fluid efflux
to permit constituent transport across the boundaries. In the event that fluid fluxes are zero
normal to the boundary, condition [3.25c] stipulates zero component flux. Caution needs to be
exercised in applying the type-2 condition on boundaries with inward fluid  velocities of
uncontaminated fluid.  If the concentration at the boundary is nonzero due to presence of
chemical near the surface, the type-2 condition will stipulate an inward mass flux density  equal to
the water phase concentration at the boundary times the phase-summed velocity, g*, which is
incorrect. In such a case, a type-3 boundary condition with zero influent concentration should be
specified to force the desired zero-flux condition. The type-2 zero gradient condition is the
default condition.

Equation [3.25d] is a type-3 boundary condition for use on boundary regions,^ which have
inward NAPL phase fluxes predicated on the assumption that contaminant enters the system
solely via inflow of nonaqueous phase liquid. Here, qon is the normal Darcy velocity of the oil
phase at the boundary  and c**Ie is the equivalent water phase concentration of species a  in the
influent oil phase defined by c**Ie =  c**1 /rao where cj" is the actual concentration of a in
the influent oil. Note that applying the condition c^Ie = 0 corresponds to imposing zero mass
flux of component a on the boundary segment. In applying the type-3 boundary condition
during periods of oil infiltration, it is important to ensure consistency in the specified oil phase
concentrations such that
                                  Cc,o
                                                                                   [3.26]

                                           23

-------
where Cso (a=l, . .  ., ns) is the oil phase concentration expressed as mass of component per
phase volume for component a.  Equivalently, in terms of mass fractions the condition may be
expressed as
                                                                                  [3.27]
where fao is the mass fraction of componenta, in the oil phase. Mass fractions and phase
concentrations are related by
                               o,papp                                         [32g]


where pp is the bulk phase density.  Simple mixing theory permits component mass fractions and
volume fractions to be related as
                                    Pa                                            [3.29]
where pa is the density of pure noninert components.  The oil phase density may be related to
phase composition by
                                I  '» /P«                                      [330]
In applying the type-3 oil infiltration boundary condition, influent oil phase concentrations,
      must b£ sPecif"ied sucn that [3.26] - [3.30] are obeyed.
3.2 Numerical Model Description

3.2.7 General solution approach

Fluid flow equations are highly coupled with each other due to interdependence of fluid
permeabilities, saturations and pressures mandating a simultaneous solution approach as
discussed in the previous chapter. The flow equations are also coupled with the transport
equations via the occurrence of interphase mass transfer terms and through the dependence of

                                           24

-------
various coefficients in the flow equations on fluid composition—e.g., phase density, viscosity,
scaling factors, etc.  Since mass transfer rates are necessarily small for the case of low solubility
organic fluids and changes in fluid properties over short time periods will accordingly also be
small, the dependence of the flow equations on transport will be very mild over short time spans
permitting computational decoupling.  Over extended periods of time,  cumulative phase mass
transfer may have significant effects on flow as dissolution and volatilization deplete the
nonaqueous phase liquid requiring some  means of updating phase transfer rates as well as
compositionally sensitive fluid properties. MOTRANS considers  fluid compositional effects on
phase density but disregards effects on other coefficients.  Interphase mass transfer rates and
phase densities are time-lagged in the flow solution and updated at the end of each time-step after
solution of the transport equations.

The transport  equations are highly  dependent on the solution of the flow equations due to the
occurrence of fluid velocity and phase saturation terms directly in the transport equations and in the
functional forms for the dispersion  coefficients.  Thus, the flow equations  must be solved
concurrently with  or prior to  .evaluating transport. Due to  the  weak back-coupling,  the  most
opportunistic approach is  to  solve  the transport-equations serially with the flow equations.
Furthermore, since the individual component transport equations are weakly coupled with each other
in the present model due to the assumption that components do not interact, the transport equations
for different components may be solved serially in arbitrary sequence.

The solution approach for solving the coupled multiphase flow and multicomponent transport
problem when equilibrium mass transfer is assumed is as follows:

1.   Solve the fluid flow equations simultaneously for the current time-step using time-lagged
    phase densities and interphase mass transfer rates.

2.   Solve the phase-summed transport equation using the previous time-step phase densities and
    interphase mass transfer rates.

3.   Back-calculate new interphase mass transfer rates and update phase densities for the current
    time step.

4.   Proceed to the next time step.

When  nonequilibrium mass transfer is considered, a modified solution approach is employed as
follows:

1.   Solve the fluid flow equations simultaneously for the current time-step using time-lagged phase
    densities and interphase mass transfer rates.

2.   Solve the phase-summed transport equation using current values of apparent partition
    coefficients, interphase mass transfer rates and phase densities.
                                           25

-------
3. Back-calculate interphase mass transfer rates, update phase densities and apparent partition
   coefficients and repeat (2) until transport solution converges.

4. Proceed to the next time step.
3.2.2 Finite element formulation

The phase-summed transport equations given by [3.15] or [3.16] are solved using an upstream-
weighted Galerkin finite element method with linear rectangular elements.  We will present the
solution for [3.15]. The formulation in terms of [3.16] is fundamentally identical. Employing
linear shape functions, Nh the concentration Cm (xt,t) can be approximated by
                                            cavj(t
                                                                                   [3.31]
where CmJ denotes the concentration at node J and e and T| are local coordinates. Detailed
information pertaining to the upstream weighting functions have been discussed by Kaluarachchi
and Parker (1989). The Galerkin finite element approximation of [3.15] over the entire domain
R can be written as
                         [ NX L (caw)  dR = 0
 [3.32]
where LfC^J is an operator for the phase-summed transport equations. A set of equations results
of the form
                                        dc
                   [P]  {Cv,,} H- [M]  {	— } - {R}
                          'Jew •
                                         d t
                                                                                   [3.33]
Using the type-3 boundary condition given by [3.25d] and disregarding cross derivative terms,
the matrices [P], [M\ and {R } are defined by
                                 d W   dN         N            dN
                         R,
                     N
                         / Wj P* Nj d R - £   < q* >    W  N  d S
                     e=l  •                 e=l
[3.34a]
                                           26

-------
                          •E
                                                                                 [3.34b]
                                                                                  [3.34c]
where TV is the total number of elements, Ns is the total number of boundary elements, W, is the
asymmetric upstream weighting function,  is the average normal oil flux along the surface S
where contaminant enters the system and < g*> is the average normal phase-summed velocity
on surface S.  If an evaporative boundary is present, the term < g*> in [3.34c] will have and
additional term (- o£*f raa) to account for the atmospheric boundary layer and the terms
containing c**Ie  in [3.34c] will vanish.

Treatment of the integrals in [3.34] employs the method of influence coefficients described in the
flow analysis and will not be repeated here. Handling of the convective term in [3.34a] using the
influence coefficient method has not been previously described and will be discussed below.
Phase-summed velocities within an element are represented by
                                                                                   r    n
                                                                                   [3.35]
where q*ig is the phase-summed velocity at the evaluation points in the element, which
may be at the nodes, at Gauss points or elsewhere as controlled by the parameter GP
(see discussion in flow analysis). Using equation [3.35] on the convective term of [3.34a]
and using linear upstream weighting shape functions (Kaluarachchi and Parker, 1989), it is
possible to show that
E  / .,
                                        3 N
                                             --> -
       N    4   /                             ,
           EX™^  I ^      xu -\     r   xcf-\
           2^,  \ — i  I H   ]g+[H   ] ff } + —-  {
      e-1   g=1  k                                                     '              [3.36]

where w and d are the width and height of a rectangular element, x and y are the coordinate axes
and matrices [//] with superscripts xg and^g refer to the contributions from standard linear shape
functions while matrices with superscripts xu and_yw refer to contributions from upstream
weighting coefficients.  Details of matrices [H] are described in detail by Kaluarachchi and
Parker (1989) and will  not be repeated here. The  surface integrals in [3.34a] and [3.34c] can be

                                           27

-------
described by
                  7r N T dS = 	  (4±3co)      if I = J
                   i  j       12

                     = J- (2 ± 3co)            if i * j                        [3.37a]
                       12
and
                            r dS =  —  ( 1 ± co)
                                                                                 [3.37b]
where Q. is the length of the side and w is the upstream weighting coefficient. The appropriate
sign of the integrals given in [3.37] will be determined by the side of the element exposed to the
boundary condition.

Time integration of [3.33]  is performed using a standard finite difference approximation. The
resulting final set of equations can be written as
 ecpi.^H-  "«:„>«. f  JfL.(e-i)  IP]   '  (ca.) •.(«)"     P-3»]
              itlJ                V     *>1                /
where k and k+l refer to previous and current time steps and 6 is a time weighting factor where
0.5 < 6 < 1 for an unconditionally stable temporal approximation. Time-step size is varied under
program control by a time incremental factor which increases or decreases time-step size
depending on the number of iterations required-for solution of the flow problem.  The time-step
size for solution of the transport problem is taken as identical to that for the flow problem. If
equilibrium mass transfer is assumed, the transport equation is linear and no iteration of the
solution is required.  For the case of nonequilibrium mass transfer, nonlinearities are handled by
a simple Picard iteration scheme with apparent partition coefficients, mass transfer rates and
phase densities updated at each iteration.  The convergence criterion for the aqueous
concentration (Caw) of species a is as follows
                     c rtl - c r
                      aw     aw
                                                                                  [3.39]
where r+1 and r refer to current and previous iterations, er is the allowable relative convergence
error which is internally set equal to 0.001 and ea is the allowable absolute convergence error for
species a. This convergence criteria must be met at all nodes.
                                           28

-------
3.2.3 Interphase mass transfer and density updating

After solution of the phase-summed transport equations, the interphase mass transfer terms Rap
must be updated. This is performed after each iteration of the transport equations in the case of
nonequilibrium mass transfer. For the equilibrium case, the transport model is linear and no
iteration is required.  Hence, for the latter case, updating is performed only once at the end of
each time step.  To compute the mass transfer rates, the concentrations of each component in the
oil, gas and solid phases are first calculated from the partition coefficients (eq. 3.13 for the
equilibrium case or 3.17 for the nonequilibrium case) and from the control  phase concentrations
(Caw for the case of eq. 3.15 or Cm for eq. 3.16). The mass transfer rate terms, Rap, are computed
by back-substitution into the component transport equations for each phase (eqs. 3.11 and 3.12)
using a finite difference approximation.  Summing over the components yields the desired ^
terms for each phase for use in the subsequent solution of the flow and phase-summed transport
equations.

For equilibrium mass transfer problems, effects of mass transfer terms on the flow and transport
equation solutions at  any given time-step will be quite small.  However, cumulative effects can
be very important since these terms will control the long term removal of NAPL associated with
dissolution and/or volatilization. The accuracy of interphase mass transfer calculations will
generally be most, troublesome during periods of highly transient NAPL flow. During such
periods, it may be preferable to suppress mass transfer updating to avoid numerical instability as
well as to reduce computational effort and improve accuracy. Therefore, MOTRANS provides
the user with the option of performing interphase mass transfer calculations at each time-step or
of not performing such calculations.  This option can be used when equilibrium mass transfer is
assumed.

With increasing simulation time, the density of each phase will change depending on the net
mass of each component leaving or entering the phase. Although density derivatives were
neglected in the equation development, on the assumption that these terms  are small over a given
time step, cumulative changes in density cannot be ignored since over long times these can
accumulate and become quite large in magnitude. To accommodate these effects, liquid phase
densities are updated at the end of each time-step as follows
                    P  = P
                    ~ w   ~ &
i-E
   a=l
                                           a=l
[3.40a]
                                  a=l
                                                                                 [3.40b]
where ns denotes the number of organic species, p° is the density of uncontaminated water, pa is
the density of pure a-component. Density ratios, pm,pro and pra, are computed as the ratio of pw,
p0and pato p°.
                                           29

-------
For the gas phase, density depends on both pressure and composition.  Pressure dependence is
treated as a nonlinear effect in the gas flow equation according to
where A is the gas compressibility taken to be 1 . 17 x 10"6 g cm'4, p° is the density of native soil
air taken to be 1 . 12 x 10"3 g cm-3, and p° is the density of contaminants in the gas phase
computed as
                                   ns
                              Pa'  = E  Cc,a
                                  «=l                                              [3.42]
where Caa is the gas phase concentration of species a and ns is the number of noninert organic
species.

In certain problems, numerical oscillations may lead to unacceptable mass balance error in the
bulk fluid phases on in the components.  A mass injection scheme has been implemented in the
code that distributes mass balance errors within the system. If this option is invoked and flow
only is being analyzed, the volumetric error for the oil phase (net change within domain minus
boundary fluxes) is computed at the end of each time step, and the error is reinjected during the
next time step in proportion to the local oil phase saturation. If flow and transport is analyzed,
the mass balance error for each component is computed at the end of each time step, taking into
account boundary fluxes of the components in water, oil and gas phases; and the error is
reinjected during the nest time step in proportion to the local oil phase saturation.
                                           30

-------
                            4. GUIDE TO USE OF MOTRANS

4. 1 Estimation of Soil and Fluid Properties

4.1.1 Soil properties

Soil properties required for analyses of multiphase flow using MOTRANS include the saturated
conductivity to water in the vertical direction, KSW , and the horizontal direction, KSW  , soil
porosity, ((), apparent irreducible water saturation, $m, maximum residual oil saturation for water
imbibition, s™*K, and the van Genuchten air-water capillary retention parameters, a and n.
Typical values for the maximum residual oil saturation, s™^x, are in the range of 0.15-0.35.
Values tend to increase for more heterogeneous porous materials and may be dependent to a
certain extent on fluid properties with higher viscosity fluids tending to have larger residual
saturations. Estimation of van Genuchten parameters from grain size distribution data (e.g.,
Mishra et al, 1989) can provide parameters which are fairly reliable, although analyses on
laboratory core samples are desirable for more accurate determinations. Estimates of hydraulic
conductivities can also be inferred from grain size distribution data, but are subject to large
uncertainty. Therefore, field or laboratory measurements to determine hydraulic conductivity is
always desirable. Typical values of capillary pressure curve parameters and vertical hydraulic
conductivities are tabulated in Table 4.1 for a variety of soil types. Horizontal conductivities are
commonly 2-5 times higher than vertical conductivities. From Table 4. 1 it is observed that ranges
of ((), Sm and n are quite narrow compared to a and Ksw. Furthermore, a useful-correlation
between the latter variables may be noted as
where A=0.5 m3 d1 ( ±50%) which provides a rough estimate of a if field data on the hydraulic
conductivity are available.

To model chemical transport, two additional porous medium parameters are required by
MOTRANS for the case of equilibrium phase transfer. These are the longitudinal and transverse
dispersivities, AL and^4r. MOTRANS assumes that these parameters may be regarded as porous
media properties which are constants for all phases and independent of phase saturation or
composition. It is known that dispersivity for nonreactive chemical transport is controlled largely
by porous media heterogeneity which typically increases with the scale of observation such that
AL typically is on the order of 1-10% of the scale of the maximum plume dimension and^4r is
several times smaller than AL.

When consideration is given to nonequilibrium phase partitioning during transport, the
nonequilibrium coefficients, koa, kwa, k^, and kws, must be specified. At the present time, little
practical guidance can be given regarding the  specification of these parameters. A few recent

                                            31

-------
laboratory studies have been reported which deal with this problem suggesting that rate constants
for homogeneous soil columns are sufficiently large to yield near-equilibrium fluid-fluid
interphase mass transfer unless velocities are unusually high. However, information on the
effects of heterogeneity at the field scale on mass partitioning is virtually nonexistent.  Since the
time scale for mass transfer between zones of differing hydraulic properties is much greater than
the time scale for pore scale mass transfer, the former may by expected to have a greater impact
on field behavior. Given site measurements of contaminant concentrations, mass transfer
coefficients may be estimated by "history matching." A better understanding of mass transfer
kinetics at the field scale awaits future study.

Table 4.1 Typical soil properties for various soil types.

                            KSW              


-------
                          Pro
                                                                                   [4.2]
where fao is the mass fraction in the oil phase of organic component a and prs is the specific
gravity of component a.

Although NAPL viscosity varies with fluid composition, MOTRANS assumes these effects are
comparatively minor and does not update viscosity for temporal phase composition changes. Bulk
phase viscosity may be determined experimentally for the fluid of concern, or an estimate may be
made from the phase composition if this is known via

                           n" = ?  f«o^ra                                       [43]
where r\ra is the relative viscosity of component a.

Capillary pressure curve scaling parameters may be estimated from surface tension and interfacial
tension data via
                                                                                  [4.4a]

                                     a
                                     °ow                                          [4.4b]
where  ow is the surface tension of pure water (i.e., 72 dynes/cm), a 0 is the surface tension of the
NAPL and  a^ is the interfacial tension between NAPL and water.

Surface and interfacial tensions of mixtures may be measured directly or estimated from the
NAPL composition via


                              O ~ £-1   010° a                                        r A  r  ~\
                                  a                                               [4  5aJ
                                                                                  r,
                                                                                  [4.5b]
                                           33

-------
where aK is the surface tension of pure a component, aao is the interfacial tension of a with
water, andfao is the mass fraction of a in the NAPL. In the absence of interfacial tension data, an
estimate of P^ may be made from the relation
                                                                                    [4.6]
Furthermore, since numerical difficulties can be induced when deviation from [4.6] occur during
the transition from an air-water to an air-oil-water system, it is recommended to employ scaling
factors whi ch force condition [4.6]. If flao and P^ are estimated from [4.4] and condition [4.6] is
not met, the condition may be imposed by the following means


                            B'   = f cor B
                              30         ao                                        [4.7a]
where

                           f
                                                                                   [4.7c]
in which B'a o and B'ow are "corrected" parameter values which may be used in lieu of B'a o  and
B'ow values which are computed directly from [4.4]. For hydrocarbons which have low solubility
in water, condition [4.6] is generally closely approximated. For gasoline, the values B'a o and B'ow
estimated from surface tension data are 3.2 and 1.45, respectively ( ± 10%).

Another empirical method for estimating scaling coefficients, based on a correlation developed
for crude oil surface tension, is
                                                                                   [4.8a]
                                       0.5
                                                                                   [4.8b]
                                       Pro
where pro is the specific gravity of the NAPL. Values of surface tensions, interfacial tensions,
specific gravities and relative viscosities for a number of common hydrocarbon compounds
which occur in groundwater are tabulated in Table 4.2.
                                           34

-------
4.1.3  Component properties

In order to model component transport, a number of additional physical and chemical properties
of the NAPL components must be known. These include the diffusion coefficients in bulk water,
oil and air (D£W , D°o and £>a°a), the oil-water, air-water and solid-water partition coefficients
(rao'  raa and ras) and the first-order decay coefficients (viaK,  vao ,  vaa and iaas).
Molecular diffusion coefficients in bulk water and air can be found tabulated for many common
industrial chemicals or they may be estimated using semi-empirical estimation methods (see
Lyman et al, 1982). Oil phase diffusion coefficients can be estimated from the water phase
diffusion coefficient by the empirical relation
                               D°
                                
-------
related the commonly used "distribution coefficient", kd, via
                                                                                   [4.12]
where pb is the porous medium bulk density. Adsorption coefficients for hydrophobic organics are
often estimated using a correlation with organic carbon content which gives
                                    foc
where/^ is the organic carbon mass fraction in the soil and ko c  is the organic carbon normalized
adsorption coefficient which may be estimated from solubility orafrom chemical structure by various
means (e.g., Lyman etaL, 1982).

The estimation of first-order decay coefficients (|iaw,  |iao, |iaa, and |ias) is more problematic due to
their nebulous physical significance. True first-order  chemical decomposition will occur due to
hydrolysis reactions which can have rate constant which can be found tabulated for many
common compounds. However, biologically mediated transformations  are generally much more
significant processes in the field and the corresponding apparent rate coefficients are dependent
on a large number of environmental variables. Since decay may be expected to occur nearly
exclusively within the aqueous phase, ignoring all terms except \im is probably justified. Existing
field and laboratory studies on the chemicals of concern under similar conditions may be used to
estimate the desired rate coefficients. However, it should be recognized that without calibration
at the site under consideration, the estimated rate coefficients will be subject to considerable
uncertainty.

A list of physical and chemical properties of selected organic compounds of environmental
significance is given in Table 4.2.
                                           36

-------
Table 4.2 PROPERTIES OF ORGANIC COMPOUNDS AT 20° C
Organic Diff. coef. Diff. coef.
compound in water in air
(m2/d) (m2/d )
Benzene
Toluene
o-xylene
Ethylbenzene
Tri chl oroethy 1 ene
Tetrachl oroethy 1 ene
1,1,1 Trichloroethane
1,1 Dichloroethane
1,3 Dichlorobenzene
Carbon Tetrachloride
1,2,4 Trichlorobenzene
9.42xlQ-5
8.21xlQ-5
6.21xlQ-5
6.21xlQ-5
7.28xlQ-5
6.56xlQ-5
T.OlxlQ-5
7.89xlQ-5
6.28xlQ-5
7.05xlQ-5
5.93xlQ-5
0.76
0.68
0.61
0.61
0.71
0.64
0.69
0.78
0.62
0.70
0.57
Relative
viscosity
(-)
0.65
0.58
0.81
0.68
0.59
0.90
1.20
0.44
1.48
0.97
2.77
Specific
gravity
(-)
0.878
0.867
0.880
0.867
1.464
1.623
1.339
1.176
1.288
1.594
1.454
Surface
tension
(dyne/cm)
27.90
29.00
31.21
29.20
29.30
31.30
25.40
24.66
37.00
27.00
45.57
Interfacial
tension
(dyne/cm)
28.88
36.10
36.06
35.48
34.50
44.40
45.00
47.20
40.00
45.00
27.13
Henry's
constan
(-)
2.40E-
2.80E-
2.20E-
3.70E-
4.20E-
3.50E-
7.70E-
2.40E-
1.10E-
9.70E-
4.30E-
                                               37

-------
4.2   Program Structure and Operation

4.2.1  Basic features and program control

4.2.1.1 Basic problem definition. MOTRANS may be used to analyze flow problems involving
two phase oil-water flow with gas at atmospheric pressure (IOIL=1, IAIR=0) or two phase air-
water including explicit consideration of gas flow (IOIL=0, IAIR=1) or three phase air-oil-water
flow (IOIL=1, IAIR=1). Problems involving flow only may be analyzed (ITRN=0) or coupled
flow and transport may be modeled (ITRN=1) for NAPLs containing up to 5 partitienable species
(NTSP). Problems may be run for plane sections (IRAD=0) or in radially symmetric sections
(TRAD=1).

4.2.1.2 I/O unit numbers used bv program. MOTRANS uses five unit numbers for data input and
output designated internally as NR, NO, NINP, and NOUT (units 5, 6, 8, and 9 respectively). NR
and NO correspond to main input and output files. NOUT is an auxiliary output unit to dump
results of the last converged time-step for restarting purposes. NINP is the unit number for the
auxiliary input data file for restarting which is identical to the file dumped by NOUT from the
previous simulation. In addition, a number of files in unformatted structure for internal storage
are used corresponding to units NWK1, NWK2 and NTRD. If the restart option is used, the
auxiliary data file must exist to define the initial conditions for the entire problem including gas
convection and transport if required. The parameter TIME for the restart problem would
normally be  adjusted to correspond to the ending time of the previous simulation.

4.2.1.3 Use of consistent units. Linear dimensions can be  in centimeters  or meters depending on
the value of IDIM and any consistent units for time may be employed. For analyses involving
transport, units of mass may be in milligrams, grams or kilograms depending on the value of
KUNIT. Units of input parameters are indicated in Appendix A and should be adhered to
religiously! Units are specified genetically as L=length, M=mass and T=time. Units  of
dimensionless variables are denoted by [-]. When consideration is given to integrated fluxes
across boundary regions, it should be kept in mind that boundary integrals from plane slice
problems will yield quantities per unit width in the direction perpendicular to the plane of the
problem, whereas radial problems are already integrated in the radial direction.
4.2.1.4 Maximum array dimensions. The program is dimensioned to accommodate problems
with up to 1500 nodes with rectangular meshes up to 50 nodes in width or height. The porous
medium may consist of up to 10 different soil types. For the main flow module, up to 100 type-1
and 100 type-2 boundary nodes or element sides may be specified for each phase described by up
to 100 schedules with up to 4 subschedules per schedule allowed. The transport module can
accommodate up to 5 non-inert species with a total of 100 type-1 and 100 type-3 boundary nodes
or element sides for each species. Transport boundaries may be described by up to 100 schedules
with as many as 4 subschedules per schedule.

4.2.1.5 Program control and termination. When data files  are first prepared or after significant
changes have been made, it is desirable to verify that the input is correct by  first reading the data
file and echoing it back. This may be done by  setting the parameter IRUN=0 which will echo

                                           38

-------
input data as well as print out information on initial conditions at all nodes. It is recommended to
always use this before running a problem. Having verified the input data to this extent, it is
suggested to first run only a few time steps to check whether convergence behavior and accuracy
(e.g., mass balance) are acceptable. Termination can be controlled for such purposes by
stipulating appropriate values of the maximum number of time-steps (MXTST) or the maximum
simulation time (TMAX). Normal program termination will occur if one of several criteria are
met:

 (1) Current time of simulation exceeds input value of TMAX
 (2) Total number of time steps exceeds input value of MXTST
 (3) Change in NAPL volume in system exceeds input value of TINF if JINF=1
 (4) Change in water volume in system exceeds input value of TINF if JINF=2
 (5) Change in gas volume in system exceeds input value of TINF if JINF=3
4.2.1.6 Convergence and time-step control. For a given phase phase p, the allowable convergence
error e;at node i for the solution of the nonlinear flow equations is given by


                           e. = h r.  e  + e
                                 PI   r    a                                      [4.14a]

and all nodes must satisfy the condition


                                    r    <:  a
                                                                                [4.14b]
where h^± and h^1 are p-phase heads at node i for the current and previous iterations, er is the
prescribed relative convergence error (RCON) and ea is the prescribed absolute convergence
error (ABSA). Termination occurs if convergence is not met in MCNT iterations. If convergence
occurs in more than MCNT/3 iterations, the simulation will proceed to the next time-step but
reduce the time-step size by a factor TFACT. If the solution converges in less than or equal to
MCNT/4 iterations, the time-step size for the subsequent time-step will be increased by a factor
of TFACT which should range between 1.01 and 1.05.

If termination occurs due to nonconvergence, the user has a number of options since many
factors influence convergence behavior. Some steps which may be taken to remedy
nonconvergence are: reduce the initial time-step size (BELT) and the timestep increment factor
(TFACT), increase the upstrealn weighting factor for the flow problem (UFW), switch off
interphase mass transfer during highly transient flow periods (KUDAT=0), apply changes in
boundary conditions gradually using "ramping" as described later, relax the relative or absolute
convergence criteria (RCON and ABSA), make the mesh finer, reduce the size of the domain, or
run the problem in stages with output from  early stages used to define later stages.
                                          39

-------
If nonequilibrium interphase mass transfer is modeled, the transport equation is solved iteratively
using a Picard scheme with convergence criteria defined by

                            &  = Cr   e  + &
                             i    ctwi   r    a                                       M ^5]

where carwi is the water phase concentration at node i for the previous iteration, er is a relative
convergence criterion which is internally set to 0.001, andea is the user-specified absolute
convergence criterion for species a in the water phase.

4.2.1.7 Upstream weighting. Upstream-weighting is a useful technique to reduce the nonlinearity
of the multiphase flow equations and obtain a more stable solution of the transport equations free
of oscillations. In selecting suitable values for upstream-weights, it is preferable to use the
minimum value necessary to obtain convergence since high values may result in numerical
dispersion at sharp fronts especially in the transport solution. Upstream-weights may range from
0 to 1 with 0 representing no upstream-weighting and 1 corresponding to full upstream-
weighting.  Values around 0.5 may be selected to start a problem and these may be adjusted
upward if oscillations in the solutions occur. A few preliminary runs may be used to evaluate the
sensitivity of the results to the upstream weights for the flow equations (UFW) and for transport
(UFT).

4.2.1.8 Automatic control of oil and water flow solutions.  MO TRANS provides two options for
controlling the oil and water flow equation solutions to reduce computational effort. In many
problems, oil saturations in the domain may reach residual levels after which they do not change
with time (except due to interphase mass transfer). If changes in oil saturation at every point in
the domain are less than a specified tolerance (TOL3) over a time-step, MOTRANS will not
solve  the oil flow equation for subsequent timesteps. That is, the system will be regarded as an
air-water system with oil in a residual occluded state. Oil saturation is updated after subsequent
time-steps to account for mass lost due to dissolution and volatilization if the transport equation
is being solved. This option is invoked by setting KFEMC=1 and specifying TOL3 > O. If
KFEMC=1 the oil equation to be solved at all times. Using TOL3 of 0.0001 to 0.001 is
satisfactory for most problems. Caution: Once the oil equation is switched off, it cannot be
reactivated unless the problem is restarted. If time-varying boundary conditions are  such that oil
may become trapped (e.g., rising water table) and later released (e.g., falling water table), the user
should be careful of using TOL3.

The second method of controlling the liquid flow equations uses the ASD method. In this
method, elements are categorized as "active" or "inactive"  at a given iteration. If the element is
active, then it is included in the global matrix assembly, otherwise it is excluded. The criteria for
changing a node from inactive to active is based on the changes in water pressure head, oil
pressure head, oil saturation, and water saturation at connected nodes from the last converged
time-step to the current iteration. The specific critenia for changing an element from inactive to
active is
                                           40

-------
                               Ahw > TOL1
                             or
                               Aho > TOL1
                             or
                               A SM > TOL2
                             or
                               A So > TOL2
where A/2W, A/z0, A5^, and A5^ are the changes in hw, h0, Sw, and S^ respectively, between the
previous converged time-step and the current iteration, TOL1 is a specified head tolerance and
TOL2 is a specified saturation tolerance. Equations for nodes on the active-inactive frontier will
have contributions only from active elements and the boundary will shift naturally as imbibition
or drainage "fronts" pass through the domain. To enable the ASD option, the switch IFEMC must
be set to 0 and TOL1 and TOL2 set to appropriate values. Values for TOL1 of 0.0001 to 0.0001
m and for TOL2 of 0.0001 to 0.001 are generally applicable for most problems.

4.2.1.9 Automatic control of gas flow solution. MO TRANS provides a number of features to
control the gas flow solution which enable gas flow to be considered while minimizing the
computational effort and improving solution robustness. The user has the option of disregarding
gas compressibility (ICOM=0) which reduces the degree of nonlinearity of the gas flow solution.
If gas pressure-differences in the domain are within a few tenths of atmospheres, or if the
transient behavior of the gas flow solution is not of concern, this option may be invoked.

In many problems of practical  interest, gas withdrawal or injection at constant pressure is
imposed which will result in a steady state flow problem after some time. To minimize
unnecessary computational effort, the user can specify a tolerance for eliminating the gas flow
equations from the flow solution when the gas pressures approach steady state. If changes  in gas
heads for all nodes in the domain are less than TMXG over a time step, gas flow will be assumed
steady and the gas equation  will  not be solved for subsequent time steps (LAIR is internally set to
0). If ITRN=1, gas velocities in the transport equation will be fixed thereafter at the steady-state
values. TMXG=0 forces continuous solution of gas flow equation.

Many gas flow problems involve domains which have large portions where the gas saturation is
zero (e.g., nodes below the water table). At such nodes; gas pressure has no physical meaning.
Solution for gas pressure at  these nodes expends fruitless effort and, more importantly, can result
in a destabilizing effect on the entire solution due to the ill-defined nature of the gas flow
equation at locations where  gas permeability is zero. If St =1 at a node, the program will
internally fix the gas pressure head to eliminate the node from the solution domain for the  gas
equation.
                                           41

-------
4.2.1.10 Mass transfer updating for transport solution. Flow and transport equations are coupled
through the interphase mass transfer terms, which are updated at the end of each time step and
time-lagged. Phase densities are updated at this time also and lagged. During periods of highly
transient flow, particularly during periods of high oil phase velocities, this decoupled approach
can lead to instabilities if errors in the back-substituted mass transfer rates grow with time. Since
the errors which result from disregarding interphase mass transfer over short time periods are
small if the component solubilities are small, it is advisable to avoid mass transfer updating
during such periods. User control over mass transfer updating is achieved through the parameter
TTRN which is the simulation time at which mass transfer updating is begun. During periods of
rapid oil infiltration, mass transfer updating may be entirely disabled by putting a large value for
TTRN. During periods of redistribution, mass transfer updating may be disabled for the first few
days. If component mass balance errors arise, it may be advisable to increase TTRN.

4.2.1.11 Mass balance correction scheme. In certain problems, numerical oscillations may lead
to unacceptable mass balance error in the bulk fluid phases or in the components. An optional
mass injection scheme has been implemented in the code that distributes mass balance errors
within the system. If this option is invoked (IMCOR=1)  and flow only is  being analyzed
(ITRN=0), the volumetric error for the oil phase (net change within domain minus boundary
fluxes) is computed at the end of each time step, and the error is reinjected during the next time
step in proportion to the local oil phase saturation. If flow and transport is analyzed (ITRN=1),
the mass balance error for each component is computed at the end of each time step, taking into
account boundary fluxes of the components in water, oil and gas phases, and the error is
reinjected during the next time step in proportion to the local oil phase saturation.
 4.2.2 Spatial discretization and mesh generation

Only rectangular elements with sides parallel to the principle axes are permitted by MOTRANS
(Figure 4.1). The domain is specified to be a rectangular region. The finite element mesh is
defined by the user by inputting z-coordinates of all nodes on a single vertical transect and x-
coordinates (or r-coordinates) on a horizontal transect.

In order to minimize the band width of the matrices which must be solved, the global node
numbering is always assumed to be in the direction of the side of the mesh with the fewest nodes,
starting in the lower left hand corner. In the event of equal numbers of nodes in the vertical and
horizontal directions, node numbering will be in the vertical direction by default. As an example,
Figure 4.2 shows a mesh with 8 rows and 9 columns. Node numbering is thus from bottom to top
as shown starting from the lower left-hand corner. Elements are numbered in the same pattern as
nodes.
Figure 4.1. Local side numbering of quadrilateral elements.
                                          42

-------
Figure 4.2.  Finite element mesh with numbers in parentheses indicating global element numbers
           and others referring to global node numbers.
It is very important to understand the global node and element numbering in order to define the
boundary conditions. For flux-type boundary conditions, it is also necessary to stipulate the local
element side number subjected to a prescribed flux. Local side numbering is always
counterclockwise from the bottom as illustrated in Figure 4.1.

When creating a mesh, several factors should be kept in mind. The mesh should be designed so
that the "aspect ratio," which is the ratio of maximum to minimum element side length, is not too
great. Generally, aspect ratios should be less than 5. The mesh should be finest near locations
where nonlinear boundary conditions will be most severe (e.g., boundaries with oil infiltration,
gas vacuum extraction, etc.). Variations in element dimensions should not occur too abruptly in
the mesh. It is recommended that changes in mesh spacing between adjacent elements are not
greater than 50%.
4.2.3  Specification of initial conditions

4.2.3.1 Initial conditions for flow. Initial conditions for the flow equations involve stipulation of
fluid pressure heads at all nodes in the domain. Initial conditions may be specified by the user in
the main data file or may be read in from an auxiliary file for restart problems. For non-restart
problems (IRES=0), initial liquid heads may be stipulated using a bilinear interpolation algorithm
with heads defined on the left and right boundaries of the domain according to the relations

                              h  = a +  b z
                              hV-c +  dz                                          [4-16]
                                                                                   [4.17]

where z is the elevation along the vertical  transect on the boundary and a, b. c and d are input
parameters. For equilibrium conditions a = zaw,  b = -\,c= pr(^ao and d = -pro where zm is the air-
water table elevation (elevation where hw=0) and zao is the airoil table elevation (elevation where
h = 0). When [4.16] and [4.17] are emploved, values  of a and c are specified on the left and right
boundaries of the domain and b and d and taken as constants for the entire domain.

Note that under equilibrium conditions with zero gauge gas pressure, the air-oil elevation in an
observation well corresponds to zao and zaw=zao-(\-pm)H0 where H0 is the oil thickness in the well.
Forcing H0=0 will produce a no-oil condition. A "no-oil" condition will exist at a node if
                               B   h  + B   h
                         ,      r ow   w   r a o  a
                         h  <.  -
                                           43
                                                                                   [4

-------
which corresponds to the condition IUFO = 1 in the data file. In this case, c and d need not be
entered and the program will assign oil heads given by equation [4.18].

MOTRANS provides the capability to restart a problem from the final time output of a previous
run.  During normal program termination, MOTRANS writes final phase pressures and historical
information needed by the fluid entrapment model to an auxiliary output file which can be read if
the restart option is invoked.  If the user specifies IRES=1, oil and water heads will be read for a
two phase problem (IOIL=1, IAIR=0).  For a three phase problem (IOIL=1,  IAIR=1), pressures of
water, oil and gas are read from the restart file if IRES=1.  If the restart file contains only oil and
water data from a previous run with (IAIR=0) but gas flow is to be considered for the restart
problem the user may use the option IRES=2 and the code will initialize gas pressures to zero. If
the problem is not a restart problem and gas flow is analyzed, MOTRANS takes the initial
conditions for the gas phase to be uniform at zero gauge pressure.

4.2.3.2 Initial conditions for transport. Initial conditions for transport are stipulated in terms of
the water phase concentrations of each  component.  Three options are provided:

    KIC=1  Species concentrations are automatically assigned initial values of zero at locations
            where oil saturation is zero, while at other locations the water phase concentration
            of each species a is assigned a default value CINITa.
    KIC=2  Initial concentrations of all species are specified for each node  in the domain in the
            main data file.
    KIC=3  Initial concentrations of all species are read in from a restart file which must have
            been dumped  from an earlier run with ITRN=1.

Note that if the problem is a restart problem for flow but the  earlier analysis did not consider
transport, KIC=3 cannot be used but KIC=1 provides a good way to continue the simulation and
take into account transport.  For problems involving NAPL addition over a relatively short period
followed by long term redistribution and migration, the initial stages of the problem may be
analyzed without considering transport with very little loss of accuracy. This can provide savings
in computational effort  and more importantly may provide a more stable and well-behaved
problem since the short term infiltration and transport problem is relatively difficult.
Consistency between the stipulated initial water phase  concentrations and the NAPL composition
requires care that the conditions given by equations [3.18]- [3.22] are fully adhered to. If
nonequilibrium transport is considered, nonequilibrium initial conditions for the restart problem
are taken into account by reading in the apparent partition coefficients for the last computed time
step in addition to the phase-summed concentration solution.
4.2.4   Specification of boundary conditions

4.2.4.1  Boundary conditions for flow model. In the flow module, type-1 boundary conditions
may be stipulated which define phase pressure heads for specified nodes, while type-2 or fluid
flux boundary conditions can be stipulated on sides of elements.  To prescribe the boundary

                                           44

-------
conditions, the user must know the global node numbers to which type-1 conditions are to be
applied and the global element numbers as well as local side numbers of boundary segments for
type-2 conditions. Global node numbering follows the pattern illustrated in Figure 4.2 with the
direction of numbering always in the direction having the smaller number of nodes. In the event
of equal numbers of nodes in the vertical and horizontal direction, the node numbering will be in
the vertical direction. Local element side numbering is illustrated in Figure 4.1.

For the flow problem, the user specifies the total number of nodes subjected to type-1 boundary
conditions for the water, oil and gas  phases (NBC(1,1), NBC(2,1) and NBC(3,1)).  For each
phase, a list of nodes  subjected to type-1 boundary conditions is given and for each node the user
indicates the number  of a schedule which contains information on the time-dependent values of
the boundary condition at the node.  A given schedule may be used for one or many nodes (if
several nodes are subject to the same boundary conditions). One list of schedules is input for the
entire flow problem without reference to the phase so a given schedule number could be used for
both water and oil flow if the occasion arose. For example, if water at node I and oil at node J
were both to be subjected to infiltration under zero pressure with the same time-dependence, then
both nodes could reference the same schedule number and the program will interpret the data in
the referenced  schedule to be water heads in one case and oil heads in the second case. Flux
boundaries are handled in a similar manner with the total number of element sides subjected to
type-2 boundary conditions for each phase defined by the variables NSS(1,1), NSS(2,1) and
NSS(3,1). For each phase, a list of elements and side numbers subjected to type-2 boundary
conditions is then given and the user indicates the number of a schedule which contains
information on the time-dependent values of the boundary condition.

Consider a simple example in which the water phase at nodes 5 and 6 of a mesh are subjected to
prescribed water pressure head defined by schedule number 1 and side 3 of element 56 is
subjected to a prescribed oil flux defined by schedule 2. Lines  17 and 18 of the main data file
would then be given as
      200       }  Line 17: NBC(I,1) for type-1 boundaries
      010       }  Line 18: NSS(I,1) for type-2 boundaries
and lines 19-20 referencing the schedule numbers and node or element numbers would be given
as
      ,      ,             i  Line 19A: list of type-1 nodes and schedules for water phase
      6      1           }

            (omitted)     }  Line 19B: list of type-1 nodes and schedules for oil phase

            (omitted)     }  Line 19C: list of type-1 nodes and schedules for gas phase

            (omitted)     }  Line 20A: list of type-2 elements, sides and schedules for water

     56      3    2        }  Line 20B: list of type-2 elements, sides and schedules for oil

            (omitted)     }  Line 20C: list of type-2 elements, sides and schedules for gas

                                          45

-------
Note that since no type-1 nodes are indicated for the oil or gas phases and there are no type-2
elements for the water or gas phases, input lines for these are omitted.

4.2.4.2. Boundary conditions for transport. Boundary conditions for the transport model are
referenced in a manner similar to that for flow, except for transport explicit input is given for
type-1 (prescribed porous media concentration) or type-3 (prescribed flux). The default
condition in the event that no explicit specification of a boundary condition is made is the type-2
condition - namely, zero normal concentration gradient which implies zero dispersive flux.  For
the type-1 condition, the boundary condition schedule will describe the time-dependent variation
in concentration at a specified node.  For the type-3 boundary condition, the code fixes the
normal flux of a given species on a boundary segment according to one of two options selected.
The type-3  "infiltration" condition takes the chemical flux equal to the product of a prescribed
influent concentration in a phase times the normal phase velocity.  The type-3 "boundary layer"
condition takes the chemical flux equal to the diffusive flux through a boundary layer of
prescribed thickness. For the infiltration condition, the schedule gives the time variation in the
influent concentration of the specified phase.  For the boundary layer condition, no schedule is
necessary since the boundary layer flux is computed entirely within the program and requires no
external specification. When employing type-3 infiltration boundary conditions, particular care
should be given to ensure consistency in the specification of influent concentrations with the
partition coefficients of the chemicals to maintain a physically meaningful system description as
discussed in the section on boundary conditions in Chapter 3. Schedules for the transport
problem are numbered independently of those for the flow problem.  That is, if there are 20 flow
schedules and 10 transport schedules, the numbering of the flow schedules will run from 1 to 10
and the transport schedules will go from  1 to 20.  The user should therefore not confuse flow
schedule 5, for example, with transport schedule 5.

4.2.4.3  Boundary condition schedules. A boundary condition schedule refers to a series of input
lines describing piece-wise linear time variations in boundary conditions for either the flow or
transport problem. Each schedule may be comprised of one or more  subschedules which
describes the  linear variation in a boundary condition over a time interval. Each subschedule
prescribes two real variables. The first variable on each line is the time of the subschedule and
the second variable is the value of the boundary condition at this time. Linear interpolation of
boundary conditions with time will be made. For example, a schedule with J subschedules will
be of the form:
                             V,
                             V2
                             V,
                     T,       VJ
The schedules are used to prescribe values of the boundary conditions for a given problem. The
V's above may represent prescribed heads for type-1 flow boundary nodes, Darcy fluxes for type-
2 flow boundary segments, porous media concentrations for type-1 transport nodes or influent
                                           46

-------
concentrations for type-3 infiltration boundary segments. The sign convention for fluid fluxes is
negative for inward flux.

In prescribing boundary condition values, it is recommended to always "ramp" values over a
small period of time rather than impose abrupt step-changes to avoid numerical oscillations that
may lead to nonconvergence or deterioration of the mass balance accuracy. For example, if a
problem involves a total simulation period of 80 days and it is desired to simulate an oil spill
infiltrating under a pressure of h0 = 0 into an initially oil-free system, then the problem involves
increasing the oil head from ho = h°ri t to h0=0 at 7=0. Note that the value of h£ri t may be
easily obtained by checking the head in the output for 7=0 when the program is run with IRUN=0.
Say h£ri c = -  0 . 67 m . To ramp the boundary conditions over a 0.1 day period the user can set
up a schedule as follows:
       3       0.0000    -0.670
               0.1000    0.00000
               10.000    0.00000

which increases the head over the 0.1 day period and then maintains a constant value. In creating
boundary condition schedules, it is important to be sure to always have the final time in the
schedule to equal or exceed the maximum simulation time plus the maximum time-step size to
ensure that the boundary condition is defined over the entire period of the simulation.
                                          47

-------
                             5. PROGRAM APPLICATIONS

Three example problems will be discussed to demonstrate the use of MOFAT. Details of the
simulations are described below. Information on computational effort required to perform the
simulations on an 80386-based microcomputer are provided in Table 5.4.

5.1 Example 1. 1-D Spill of Two-Component Hydrocarbon

The first example problem involves a spill of a mixture of toluene and o-xylene at the soil
surface followed by gradual leaching by rainfall.  The domain of the problem is a 200cm vertical
section with a water table at a depth of 150 cm described by a single column of forty 5 x 5 cm
elements using 82 nodes.  Properties of the porous medium are given in Table 5.1. Note that
since the problem is one dimensional, K   alone controls the flow.  However, since the solution
        r                         '  swz                               '
employs two dimensional elements Ksv  is still  employed.  To minimize the possibility of
oscillations in 1-D vertical problems, it is advisable to input Ksv > Ksv . Here, we use Ksw  = 2
Ksv The nonaqueous liquid entering the system is assumed to be composed of 50% mass
fractions of toluene and o-xylene. The oil specific gravity, relative viscosity and scaling factors
of the mixture estimated by [4.2] - [4.7]  and the data in Table 4.2 are given in Table 5.1.
Component values  of air, oil and water phase diffusion coefficients, and oil-water and air-water
partition coefficients  estimated from [4.9] - [4.11] and the data in Table 4.2 are also summarized
in Table 5.1. Solid-water partition coefficients and decay coefficients were assumed to be zero
for both components. The problem is simulated in three stages involving one initial run and two
restart runs as follows:

            Stage  1:  Constant bead oil infiltration
            Stage  2:  Redistribution with transport for 25 days
            Stage  3:  Constant flux water infiltration for 100 days

In Stage 1, oil infiltration is considered into soil which is initially oil-free and in equilibrium with
a water table at a depth of 150 cm.  Oil is added under a constant head of h=Q cm at the top two
nodes under the JINF=2 option until the cumulative oil infiltration is 21 cm3 (=4.05 cm3cm'2).
Other boundaries for  oil are no flow. The bottom nodes are maintained at a constant water head
of hw=50 cm corresponding to the initial conditions throughout Stage 1 and other boundaries are
no flow for water.  Transport is not considered  during Stage 1 (ITRN=0) since the total duration
of the first stage is only 0.0074 d At the end of infiltration, the oil front is at a depth of
approximately 20 cm  (Figure 5.1).

Stage 2  involves a 25 day period during which oil redistribution occurs with no flow conditions
at the upper boundary. Initial conditions for flow are read from the restart file. Transport is
considered during Stage 2 with initial conditions computed internally using the KIC=1 option.
Initial oil phase concentrations of toluene and o-xylene computed from [3.20] corresponding to
50% mass fractions are 436.Img cm'3. Equilibrium water phase concentrations (CINIT) for
toluene  and xylene computed from [3.13a] are 0.2583 and 0.0762 mg cm'3,  respectively.  Stage 2
is carried out for a period of 25 days with no flow boundary conditions for  oil and with boundary

                                           48

-------
conditions for water as in Stage 1. Type-3 boundary conditions for transport are employed with
the influent concentration, c"l\ equal to zero corresponding to zero component mass flux on all
boundaries.  The water and total liquid saturation profiles at the end of Stage 2 are shown in
Figure 5.1. Note that a residual oil saturation of about 7% occurs in the upper unsaturated zone
due to a limitation on drainage under gravity when the oil permeability is low.  Oil permeability
increases near the water table as the water saturation increases resulting in faster oil drainage near
the water table.
Table 5.1 Input parameters for Example 1.
Soil properties:
Ksw = 800 cm dl
K " = 400 cm dl
sw
$ ' =0.4
Sm =0.05
Sor =0.2
a = 0.05 cm1
n =2.5
AL = 2.0 cw
AT = 0.2 cm

Bulk fluid properties:
P_2.1 (3^=1.83
rjro = 0.695 pro = 0.873

Component properties:
o-Xylene Toluene
Da°w = 0.620 0.821 cm2dl
D°0 = 0.70 0.987 cm2dl
D° = 6099. 6765. cm2dl
a a
rKO = 5729. 1683.
Taa = 0.22 0.28
pa = 880. 862. mg-cW3
 Figure 5.1 Water and total liquid saturation vs depth at end of Stages 1, 2 and 3 for Example 1.
At the beginning of Stage 3, water infiltration is initiated at the upper surface at a constant flux of
qw =10 cm dl via a type-2 flow boundary condition.  The water head was maintained as in the
previous stages at 50 cm at the bottom of the column. Zero oil flux was assumed on all
boundaxies. Initial conditions for flow and transport were read from a restart file from the end of
Stage 2 (IRES=1 and KIC=1). Assuming infiltration of clean water, a type-3 transport boundary
condition with cfv = 0 was imposed at the upper surface. At the lower surface and elsewhere a
zero concentration gradient (type-2 condition) was employed to allow convective mass loss in the
water phase from the bottom of the column.

Several different simulations of Stage 3 were run using both equilibrium and nonequilibrium
interphase mass transfer analyses. To provide greater mass balance stability for the
nonequilibrium analyses, the oil-based  phase-summed transport formulation was employed.  In
                                           49

-------
addition to an equilibrium simulation, analyses were performed using rate constants for all phase
pairs of 1.0,  0.01 and 0.005 dl.  Concentrations of xylene and toluene at the lower boundary are
shown as a function of time for the various runs in Figure 5.2. Due to its higher solubility,
toluene is lost from the system more rapidly than xylene leading to an increase in the mass
fraction of xylene in the oil phase, and hence an increase in aqueous phase concentration at short
times, followed by a reduction as extraction in the water phase continues.  As the phase transfer
rate constants increase to 1.0 dl, the nonequilibrium model is seen to give results identical to the
equilibrium model. This degeneration to the equilibrium case provides verification of the
nonequilibrium model formulation.  As the rate coefficients decrease, the exit concentrations,
and rate of leaching of contaminants from the oil phase, both diminish as expected.
Figure 5.2   Aqueous xylene and toluene concentrations at outflow boundary for Example 1 for
            equilibrium and nonequilibrium analyres. Rate coefficients in dl.
5.2  Example 2. 2-D Planar Spill with Sloping Water Table

This problem involves a light hydrocarbon spill in a two dimensional Cartesian section through
saturated and unsaturated zones with a water table gradient.  The problem is simulated in two
stages involving one initial run and one restart run as follows:

       Stage 1: Constant head oil infiltration on strip
       Stage 2: Redistribution with transport for 25 days

The domain is 8 m in the vertical and 11 m in the horizontal with a water table initially located 4
m from the bottom on the left boundary and 3.5 m on the right. The lower surface is an aquitard
and the top is the soil surface. An oil leak occurs on a strip 2 m in width centered 5 m from the
left boundary along the top.  Spatial discretization is achieved with 88 elements of 1 x 1 m
spacing for a total of 108 nodes.  Initially the water table is assumed to vary linearly from the left
to right boundaries and water pressures are hydrostatic with no oil in the system.  Boundary
conditions for the water phase in Stage 1 are no flow on the top and bottom boundaries and
prescribed head on the side corresponding to the initial conditions allowing water flow in the
aquifer from the left to the right. Boundary conditions for oil are no flow except on the strip
source. (Note: this presumes the physical domain is large enough that free oil never reaches the
boundaries.) At the beginning of Stage 1, oil infiltration commences on the strip source under a
constant head of h0 = - 0.1 m and is continued until the total accumulation of oil is 1.0 m3 (per m
in the third dimension).  The organic liquid is assumed to be 89.5% inert oil and 10.5% benzene.
Soil properties and bulk fluid properties are given in Table 5.2. Transport is not simulated during
Stage 1 which lasts for only 4.17 d.

During Stage 2 oil redistribution and benzene transport is simulated subject to the same boundary
conditions for the water phase as in Stage 1 and with no flow for oil on all boundaries. A zero
concentration gradient was imposed on the downstream boundary to allow mass loss with the

                                           50

-------
water phase by convection.  The zero gradient condition was employed on the upper and lower
surface which have zero fluid flux resulting in the equivalent of a zero mass flux condition for
transport.  On the upstream boundary, a type-3 condition with zero influent concentration was
imposed.  Stage 2 was continued for a period of 25 days during which time the solution indicated
a mass loss of benzene due to dissolution and transport of about 2%.  Contours of oil saturation
and water and gas phase concentrations at the end of Stage 2 are shown in Figures 5.3-5.5.

Table 5.2 Input parameters for Example 2.
Soil proper ties:                     Bulk fluid properties:
K
4> '
sm
sor
a
n
AL
AT
= 5.0cmdl
' =0.35
= 0.05
= 0.2
= 5.0m'1
= 2.8
= 0.2m
= 0.04m
rjro = 2.0 pro = 0.832
Component properties (benzene)
D°w = 0.94xlQ-4 m2dl
Dag = 1.13x10 m d
D° = 0.763 m2dl
a a
F = 493
-*- KO t^J.
raa = 0.24
pK = 0.87xl06 gm3
Figure 5.3 Oil saturation contours at the end of Stage 2 for Example 2.


Figure 5.4 Water phase concentration contours (gm3) at the end of Stage 2 for Example 2.


Figure 5.5 Gas phase concentration contours (gm'3) at the end of Stage 2 for Example 2.


5.3 Example 3. 2-D Radial Dense Solvent Spill with Vacuum Extraction

This problem involves a spill of tetrachloroethylene (PER) in a radial domain and remediation
using vacuum extraction. The simulation is performed in three stages with two restarts as
follows:

       Stage 1: Oil infiltration event
       Stage 2: Redistribution and transport under natural gradients
       Stage 3: Remediation using gas vacuum extraction
                                           51

-------
The first stage of the problem involves infiltration PER on a circular area with a radius of 2.0 m
at the soil surface.  A total of 3 m'3 of NAPL is assumed to infiltrate under a water-equivalent oil
head of-0.25 m.  A water table occurs at a depth of 3.2 m and an impermeable layer occurs 5 m
below the surface.  The soil is uniform and has properties given in Table 5.3. Fluid properties for
PER are also given in the table.  The problem is analyzed as a 2-D radial section with an inner
radius of 0.1 m (in order to facilitate subsequent analysis of pumping from a well of the same
radius) and an outer radius of 8 m. A mesh with 12 nodes in the vertical direction and 16 nodes in
the horizontal direction is employed.

Initial conditions for the water phase are assumed to correspond to equilibrium with the water
table with no oil. Boundary conditions for the water phase involve a type-1 condition on the
right side below the water table corresponding to vertical hydrostatic conditions (i.e., same as the
initial conditions).  All other boundaries are no flow for the water phase.  Boundary conditions
for the oil phase  are a type-1 condition with h0 = - 0.25 m on the 5 nodes on the upper surface
between r = 0 and r = 2.0 m. All  other nodes are no flow for oil. Gas flow was not considered
during infiltration (NPH=2). Termination of Stage 1 occurred at t = 6.48 d after infiltration of
3.128m3 of oil.

Stage 2 of the problem involves  a continuation of the Stage 1 using the restart option (IRES=1).
Redistribution of the PER is permitted for a period of 25 d under no flow boundary conditions
for the oil phase  on all boundaries.  Boundary conditions for the water phase are maintained as in
the infiltration stage. Gas phase flow is again disregarded. Transport is considered with the
KIC=1 option and the initial aqueous concentration of PER at nodes with oil phase set to the
solubility of 150 g m'3.  Boundary conditions for transport are type-2 (zero dispersive flux) on all
boundaries
which is equivalent to zero flux since there is essentially no flow on the boundaries during
Stage 2. By the end of the redistribution period, the NAPL plume has reached the lower aquifer
boundary.  However, a substantial volume of the spill is still retained in the zone above the water
table as "residual saturation" retained by capillary forces. The oil saturation distribution at the
end of Stage 2 is shown in Figure 5.6.

In Stage 3,  remediation is simulated with vacuum pumping in the unsaturated zone. A vacuum
well with a 0.1 m radius is assumed to be placed at the left boundary screened over a 2.5 m
interval from a depth of 1.25 m to 3.25 m and regulated at a pressure head ofha = - \.5rn. The
gas boundary at the well is treated as a type-1 boundary condition. Gas inflow is permitted along
the upper surface from r =3.5 m to the outer perimeter under a type-1 condition for gas with  a
constant pressure ofha=0. The inner 3.5m of the surface, which is assumed covered, and other
boundaries  are treated as no flow for the gas phase. Water head is prescribed on the right
boundary below the water table as in the earlier simulations and all other boundaries are treated
as no flow for water. No flow conditions are imposed for the oil phase on all boundaries.
Boundary conditions for transport are type-3 with zero influent concentration on the top and right
side boundaries and zero dispersive flux elsewhere.  Gas pumping results in water table
upwelling as shown in Figure 5.7.  The gas flow rate stabilizes at 319 m3 d1 after 2 days and the
PER recovery rate reaches a corresponding steady rate of 16.11 kg d1. PER mass in the soil  vs

                                           52

-------
venting time (Figure 5.8) shows a small increase ( ~ 1.5%) during the first day due to numerical
error before the solution stabilizes and mass removal rate becomes nearly constant.
Table 5.3 Input parameters for Example 3.

Soil proper ties:                     Bulk fluid properties:
     K  = S.Ocmd1                  Pfl0-2.3          B  = 1.77
       s w                              r ao -             r ow
     Ksw'=2.0cmdl                  r|,0 = 0.9         pro=1.62
     cj)  =0.4                     Component properties (benzene)
     Sm  =0.1                        D°w = 0.85xlQ-4 m2dl
     Sor  =0.2                        D;O = 0.95xlO-4 m2dl
     a  =3. Om-1                     D°a = 0.65     wV1
     n   =2.8                        rao  = 10,800.
     AL  =0.lm                      Taa  = 0.35
     AT  =0.02m                     pa  = 1.62xl06 gm3
Figure 5.6 Oil saturation contours at the end of Stage 2 for Example 3.


Figure 5.7 Water saturation contours at the end of Stage 3 for Example 3.


Figure 5.8 PER mass in system vs time during Stage 3 of Example 3.
                                          53

-------
Table 5.4 Computaional effort required for example problems.
Problem
Example 1
Stage 1
Stage 2
Stage 3
Example 2
Stage 1
Stage 2
Example 3
Stage 1
Stage 2
Stage 3
CPU time
(min)*

0.5
5.8
13.8

12.3
21.3

118.
139.
281.
Total number
of time steps

13
358
858

211
327

203
327
260
Avg. iterations
per time step

5.5
2.2
1.7

1.0
1.0

3.5
1.3
3.0
        CPU time for 33 Mhz 80386 computer with Weitek 3167 coprocessor running under
        DOS compiled using NDP FORTRAN
                                        54

-------
                                   6.  REFERENCES
Bear, J. 1972. Dynamics of Fluids in Porous Media. American Elsevier, New York.

Huyakorn, P. S. and K. Nilkuha. 1978.  Solution of transient transport equation using an
upstream weighted finite element scheme.  Apple. Math. Modeling. 3:7-17.

Kaluarachchi, J. J. and J. C. Parker. 1989.  An efficient finite element method for modeling
multiphase flow in porous media. Water Resour. Res. 25:43-54.

Kaluarachchi, J. J. and J. C. Parker. 1990.  Modeling multicomponent organic chemical transport
in three fluid phase porous media.  J. Contam. Hydrol.  In press.

Kaluarachchi, J. J. and J. C. Parker. 1990.  Effects of nonwetting liquid entrapment on
multiphase flow in porous media. Submitted to Trans, in Porous Media.

Land, C. S. 1968. Calculation of imbibition relative permeability for two and three phase flow from
rock properties. Trans. Am. Inst. Min. Metall. Pet. Eng. 243, 149-156.

Lyman, W. J., W. F. Reehl and D. H. Rosenblatt. 1982. Handbook of Chemical Property
Estimation Methods.  McGraw-Hill, New York.  960 p.

Millington, R. J. and J. P. Quirk. 1959.  Gas diffusion in porous media. Science.  130:100-102.

Mishra, S., J. C. Parker and N. Singhal. 1989. Estimation of soil hydraulic properties and their
uncertainty from particle size distribution data. J. Hydrology. 108:1-18.

Parker, J. C., R. J. Lenhard and T. Kuppusamy. 1987. A parametric model for constitutive
properties governing multiphase flow in porous media. Water Resour. Res. 23:618-624.

Parker, J. C.  1989.  Multiphase flow and transport in porous media. Reviews of Geophysics.
27:311-328.

Sleep, B. E.  and J. F.  Sykes. 1989.  Modeling the transport of volatile organics in varibly
saturated media.  Water. Resour. Res. 25:81-92.

van Genuchten, M. Th. 1980.  A closed-form equation for predicting the hydraulic conductivity
of unsaturated soils. Soil Sci.  Soc. Am. 44:892-898.
                                           55

-------
                APPENDIX A. DESCRIPTION OF MAIN INPUT FILE
1.    Input Data for Flow Module
Line     Variable
 Format   Description
         TITLE
72A1      Description or title of the problem
2A
          IRUN
          IRES
15
15
          IRAD
         IOIL
          IAIR
          ITRN
          IDIM
15
15
15
15
15
Index for execution control;
 = 0 reads and echoes input data
 = 1 executes the problem

Restart option:
= 0 not a restart problem
= 1 restart problem for number of phases
    same in previous and current problems
= 2 restart problem for IAIR=0 in previous
    problem and IAIR=1 in current problem
    while IOIL=1 in both (see 4.2.3.1).

Index to describe 2D flow problem:
 = 0 for plane section in x-z coordinates
 = 1 for radial section in r-z coordinates

 Index to  describe if oil phase is active,
   = 0 oil phase equation not solved
   = 1 oil phase equation is solved (see 4.2.1.1)

Index to describe if phase is active:
   = 0 gas phase equation not solved
   = 1 gas phase equation is solved(see 4.2.1.1)

Index to include transport module:
   = 0 flow only
   = 1 flow and transport

Index to specify the units of linear dimensions:
   = 1 for centimeters
   = 2 for meters (see 4.2.1.3)
                                          56

-------
Note:

2B
Note:

2C
  IFEMC         15          Index to activate ASD algorithm:
                               = 0 ASD algorithm active
                               = 1 ASD not active (see 2.2.4)

  KFEMC         15         Index to inactivate oil phase equation
                               = 0 oil phase equation active
                               = 1 oil equation inactive if A 5^ < TOL3
                                (see 4.2.1.8)

Line 2B should be entered only if IFEMC = 0

 TOL1           E10.4       ASD tolerance for liquid phase pressures

 TOL2           E10.4       ASD tolerance for liquid phase saturations
                                  (see 2.2.4)

 Line 2C should be entered only if IFEMC=1 and KFEMC=1
 TOL3
         MXTST
         MCNT
E10.4     Tolerance for oil saturation for eq. elimination
            (see 4.2.1.8)
                  15
                  15
NMAT
ICN
IROW
KCOL
KEQUL
15
15
15
15
15
         IMCOR
                  15
          Maximum number of time steps allowed
            (see 4.2.1.5)

          Maximum number of iterations allowed

          Number of material types (max=10)

          Number of flow boundary condition schedules
          (max=100)

          Number of rows in the mesh

          Number of columns in the mesh

          Index to describe interphase partition:
          = 0 local equilibrium between phases assumed
          = 1 nonequilibrium interphase mass partitioning

          Index to activate mass  injection scheme:
             = 0 mass injection scheme inactive
             = 1 mass injection scheme active
                (see Section 4.2.1.11)
                                         57

-------
                    ISOLVE
                    ITPNT1
                 15
                 15
           Index to control matrix solution method:
             = 1  Solver A - more efficient but less robust
             = 2  Solver B - less efficient but more robust

           Index for intermediate printout
            = 0 results are printed at all nodes
            = 1 results printed only at selected nodes
           Note:  The maximum number of nodes (IROWx KCOL) must not exceed 1500 and the maximum
           difference between node numbers of adjacent nodes must not exceed 50.
           4A
            4B
           6A
ZCO(I)
I=1,IROW

XCO(I)
1=1,KCOL

JINF
TINF



TIME


DELT

TMAX

DTMX

TFACT

TPNT
7E10.4
7E10.4
                                     15
15



E10.4


E10.4

E10.4

E10.4

E10.4

E10.4
z-coordinates of the nodal points on the z-axis
starting from the lower left corner [L]

X-coordinates of the nodal points on the x-axis
 starting from the lower left corner [L]

Index for terminating the current run based
on total change in phase volume in system:
= 0 option inactive
= 1 Water volume prescribed
= 2 oil volume prescribed
= 3 gas volume prescribed

Cumulative fluid intake (see 4.2.1.5) [L3]
program halts when phase accumulation>TINF
(needed only if JINF >0)

Starting time of the simulation. Usually zero
unless a restart problem [T]

Starting time increment [T]

Maximum simulation time [T], (see 4.2.1.5)

Maximum allowable time increment [T]

Incremental factor for time-step (1.02-1.05)

Time interval for printout of results [T)
Note:      Lines 6B and 6C should be entered only iflTPNTl = 1

           6B        ITPNT2        15          Number of nodes where results are to be printed

                                                     58

-------
6C
ITPNT(I)        1015
 I=1TTPNT2
 Nodes at which results are to be printed
           VISOW
           RHOW
           BAG
           BOW
           TH
           ABSA
           RCON
           UWF
                  El0.4      Ratio of oil to water viscosity, r\m [-]

                El 0.4        Ratio of oil to water phase density, pro [-]
                             (not used ifITRN=l)

                El0.4        Air-oil phase scaling parameter, Pao [-]

                E10.4        Oil-water phase scaling parameter, p^ [-]
                             (see 4.1.2)

                E10.4        Time weighting factor, 0.50<9< 1.0
                             6=1 refers to fully implicit scheme

                E10.4        Absolute convergence limit for all fluid heads [L]
                             (see 2.2.3)

                El0.4       Relative convergence limit for all phases [-]
                             (see 2.2.3)

                E10.4       Upstream weighting coefficient for flow
                            (see 2.2.1)
Note: Lines 9 to 1 ID are entered only if IRES=0
          IUFW
          IUFO
                15
                15
Index describing initial condition for water:
   = 0 water pressure from bilinear interpolation
   = 1 nonuniform water pressure distribution

Index describing initial condition for oil
   = 0 oil pressure from bilinear interpolation
   = 1 no oil condition
   = 2 nonuniform oil pressure distribution
Note:   Line 10A should be entered only if IUFW = 0

10A       WPI(l)         E10.4       Parameter a in eq. [4.16] for left boundary

           WPI(2)         El0.4       Parameter a in eq. [4.16] for right boundary
                                                 59

-------
           CWC           E10.4       Parameter bin equation [4.161


Note: Line 10B should be entered only iflUFO = 0

10B       OPI(l)            E10.4     Parameter c in eq. [4.17] for left boundary

           OPI(2)            E10.4     Paxameter c in eq. 14.17] for right boundary

           COC            E10.4       Parameter Jin equation [4.17]

Note:  Line 11C should be entered only iflUFW = 1

  11C      CHAF(1,1)      7E10.4      Initial water heads entered sequentially for all
            1= 1 ,NPT,NPH              nodes from node 1 to NXY=IROW x KCOL [L]

Note: Line 11D should be entered only ifIUFO=2

  1 ID      CHAF(I,1)      7E10.4      Initial oil heads entered sequentially for all
            I=2,NP, NPH               nodes from node 1 to NXY=IROWxKCOL [L]

Note: Line 12 should be repeatedNMAT times

  12        PROP(5,I)     E10.4        Km in the x-direction for soil I [L T1]

            PROP(6,I)     E10.4        Km in the z-direction for soil I [L T1]

            PROP(4,1)     E10.4        Total porosity, <|>, for soil I [-]

            PROP(2,I)     E10.4        Residual Water saturation, Sm, for soil I [-]

            PROP(8,I)     E10.4        Max. residual oil saturation, Son for soil I [-]

            PROP(1,I)      E10.4        Parameter a of VG model for soil I [L4]

            PROP(3,I)     E10.4        Parameter n of VG model, for soil I [-]


Note:   Line 13 should be entered only if NMAT > 1

13          IMAT         15            Index for assiong material types to nodes:
                                         = 1 material type read for blocks of nodes
                                         = 2 material types read for each node
                                               60

-------
 Note:  Line 14 should be entered only iflMAT = 1

 14A        KBL           15            Number of blocks of nodes


 Note:  Line 14B should be read KBL times
 (important: all nodes must be included in blocks!)

 14B        Kl           15             Lowest node number in block J

            K2          15             Highest node number in block J

            K3          15             Material type for block J

 Note:   Line 15 should be read only iflMAT = 2

 15     IPROP(I)     15                Material type for node I (enter sequentially for
        1= 1 ,NXY                       all nodes from node  1 to NXY=IROWxKCOL)

 Note:  Line 16 should be read only iflOIL=1 and I AIR=1

16        ICOM           15          Index describing gas compressibility:
                                         = 0 compressible
                                         = 1 incompressible  (see 4.2.1.9)

          TMXG           E10.4      Tolerance for deciding if gas flow is at steady-
                                      state. If changes in gas heads for all nodes in the
                                      domain are less than TMXG over a time-step,
                                      gas flow will be assumed steady and the gas
                                      equation will not be solved for subsequent time-
                                      steps (IAIR is set to 0). If ITRN=1, gas
                                      velocities in the transport equation will be fixed
                                      thereafter at the steady-state values. TMXG=0
                                      forces continuous solution of gas flow equation.
                                      (see 4.2.1.9)

 Note:  Total numbers oftype-l and type-2 nodes for each phase must not exceed 100.

 17        NBC(1,1)     15           Total number  of nodes with type-1 BC for water

           NBC(2,1)     15           Total number  of nodes with type-1 BC for oil

           NBC(3,1)     15           Total number  of nodes with type-1 BC for gas
                                      (see 4.2.4.1)

                                                61

-------
 18        NSS(1,1)      15          Total no. of elements with type-2 BC for water

           NSS(2,1)      15          Total no. of elements with type-2 BC for oil

           NSS(3,1)      15          Total no. of elements with type-2 BC for gas
                                     (see 4.2.4.1)


 Note: Line 19A is entered only ifNBC(l,l) > 0 and is repeated NBC(1,1) times

 19A      NNB1(2,L,1)     15        Number of the node with type-1 BC for water

           NNB1(1,L,1)     15         Schedule number for this node
           L=1,NBC(1,1)

 Note: Line 19B is entered only ifNBC(2,1) > 0 and is repeated NBC(2,1) times

 19B      NNB2(2,L,1)      15        Number of the node with type-1 BC for oil

           NNB2(1,L, 1)     15         Schedule number for this node
           L=1,NBC(2,1)

 Note: Line 19C is entered only ifNBC(3,1) > 0 and is repeatedNBC(3,1) times

 19C      NNB3(2,L,1)      15         Node number of the node with type-1 BC for gas

           NNB3(1,L,1)     15          Schedule number for this node
           L=1,NBC(3,I)

 Note: Line 20A is entered only ifNSS(l,l) > 0 and is repeatedNSS( 1,1)  times

  20A     NNSW(2,L,1)    15         Number of element with type-2 BC for water

           NNSW(3,L,1)    15          Side of the element exposed to flux BC

           NNSW(1,L,1)    15          Schedule number for this element side
           L=1,NSS(1,1)

 Note: Line 2OB is entered only ifNSS(2,1) > 0 and is repeatedNSS(2,1)  times

20B       NNSO(2,L,1)     15         Number of the element with type-2 BC for oil

          NNSO(3,L,1)     15         Side of the element exposed to flux BC
                                               62

-------
          NNSO(1,L,1)    15         Schedule number for this element side
          L=1,NSS(2,1)

 Note:  Line 20C is entered only ifNSS(3,1) > 0 and is repeatedNSS(3,1) times

 20C      NNSA(2,L,1)    15         Number of the element with type-2 BC for gas

          NNSA(3,L,1)    15         Side of the element exposed to flux BC

          NNSA(1,L,1)    15         Schedule number for this element side
          L=1,NSS(3,1)

 Note:  Lines 21 and 22 are repeated as a block ICN times to describe the ICNflow boundary condition
 schedules. For each schedule, line 21 should be given once and line 22 is repeated ISUB(L)-1 times. If
 there is only one subschedule for a given schedule (ISUB(L) = 1 then line 22 is omitted for that
 schedule.

21        ISUB(L)          15         Number of subschedules in schedule L (max=4)

          CN(1,K,L)        E10.4     Time for subschedule K of schedule L

          CN(3,K,L)        E10.4     Value for subschedule K of schedule L

 22       CN(1,K,L)        E10.4     Time for subschedule K of schedule L

          CN(3 K,L)        E10.4     Value for subschedule K of schedule L
          K=1,ISUB(L)               (see 4.2.4.3)
          L=1,ICN
  2. Input Data for Transport Module (enter if ITRN=1)

 23         KUNIT         15          Index to specify mass units:
                                          = 1  milligrams
                                          = 2 grams
                                          = 3  kilograms (see 4.2.1.3)

            NTSP          15          Number of partitienable organic species to be
                                        simulated (max= 5)

           NCON(l)        15          Total number of type-1 nodes for transport
                                        (max=100), (see 4:2.4-2)
                                               63

-------
 24
 Note:

25
  NEBC(l)         15          Total number of type-3 elements for transport
                               (max=100)

  KCN             15          Total number of schedules describing transport
                               boundary conditions (max=100)

  IEVAP          15           Index to specify presence evaporative boundary
                                = 0 evaporative boundary not present
                                = 1 evaporative boundary is present

  TRANS(l)       E10.4         Longitudinal dispersivity [L]

  TRANS(2)       E10.4         Transverse dispersivity, (see 3.1.1) [L]

  UWT            El0.4         Upstream weighting factor for transport [-]

  TTRN           E10.4         Time at which interphase mass transfer update
                                starts, (see 4.2.1.10) [T]

Line 25 should be repeated NTSP times

  GAM(a,l)       E10.4
 Equilibrium partition coefficient between oil and
 water for species, (see 3.1.2) a[-]
                            E10.4
                            E10.4
 Note:
  GAM(a,2)


  GAM(a,3)


  DCAY(a,l)    E10.4


  DCAY(a,2)    E10.4


  DCAY(a,3)    E10.4


  DCAY(a,4)


Line26A should be repeated NTSP times
 Equilibrium partition coefficient between air and
 water for species a [-]

 Equilibrium partition coefficient between solid
 and water for species of a [-]

First order decay coefficient for species a in
water phase, (see 3.1.1) [T"1]

First order decay coefficient for species a in
oil phase [T"1]

First order decay coefficient for species a in
air phase [T"1]
                            E10.4     First order decay coefficient for species a on
                            solid phase [T"1]
                                                 64

-------
26A      DIFW(a)
                E10.4
DIFO(a)
DIFA(a)
DLIQ(a)
E10.4
E10.4
E10.4
           Diffusion coefficient of species a in water [L2T4]
           (see 3.1.1)

           Diffusion coefficient of species a in oil [L2T4]

           Diffusion coefficient of species a in air [L2!"1]

           Liquid density of pure species a [M L"3]
Note:  Line 26B should be repeatedNTSP times ifKEQUL=l

26B       EPART(a,l)     E10.4
 27
EPART(a,2)

EPART(a,3)

EPART(a,4)

ABSC(a)



 KIC
            CINIT(K)
            K=1,NTSP
                            Water-gas mass transfer coeffident [T1]
                            (see 3.1.3)
                            E10.4     Oil-water mass transfer coefficient [T"1]

                            E10.4     Oil-gas mass transfer coefficient [T1]

                                       Water-solid mass transfer coefficient [T"1]
E10.4

E10.4



 15
Absolute convergence criteria for transport
eqn in terms of water concentration [M L"3]
(see 4.2.1.6)

Option for initial conditions for transport:
= 0 assign zero concentration were S0 = 0
    and specified constant values where S0* 0
= 1 read initial concentrations from restart file
= 2 read initial concentrations for all
    nodes from main data file (see 4.2.3.2)
                   5E10.4     Initial water phase concentration for species K
                              where S0* 0 (needed if KIC=0 and IUFO* 1)
Note:  Line 28 is entered only if KIC=2 and is repeated NTSP times.

28
 CON(I,cc,l)
 I=1,NXY
 7E10.4     Nonuniform initial water phase concentration of
             species a input serially for all nodes [M L"3]
                                                65

-------
 Note:  Line 29 should be entered only ifNCON(l) > 0 and is repeated NCON(l) time.
 OmitJSCH(J,l, a) entries on each line for values of a > NTSP.

 29         NSPB(J, 1)       15          Number of the node with type-1 BC

            JSCH(J,1,1)      15          Schedule number for this node for species 1

            JSCH(J,1,2)      15          Schedule number for this node for species 2

            JSCH(J,1,3)      15          Schedule number for this node for species 3

            JSCH(J,1,4)      15          Schedule number for this node for species 4

            JSCH(J,1,5)      15          Schedule number for this node for species 5
            J=1,NCON(1)                (see 4.2.4.2)


 Note:  Line 30 should be entered only ifNEBC(l) > 0 and is repeated NEBC(l) times.
 Omit KSCH(J, 1, a) entries on each line for values of a > NTSP.

 30        IDBC(J,1,1)        15           Number of element with type-3 BC (see 4.2.4.2)

           IDBC(J, 1,2)        15           Side number of the element with type-3 BC

           KSCH(J,1,1)       15           Schedule number for species 1

           KSCH(J, 1,2)       15           Schedule number for species 2

           KSCH(J,1,3)       15           Schedule number for species 3

           KSCH(J, 1,4)       15           Schedule number for species 4

           KSCH(J,1,5)       15           Schedule number for species 5
           J=1,NEBC(1)

 Note:  Lines 31 and 32 are repeated as a block KCN times to describe the KCN transport boundary
 condition schedules. For each schedule, line 31 should be given and line 32 is repeated KSUB(L)-1
 times.  If there is only one subschedule for a given schedule (KSUB(L) = 1) then line 32 will be omitted
 for that schedule (see 4.2.4.3).

31        KSUB(L)         15         Number of subschedules in schedule L (max=4)

          ACN(1,K,L)      E10.4      Time for subschedule K of schedule L
                                               66

-------
         ACN(3,K,L)      E10.4     Value for subschedule K of schedule L
32       ACN(1,K,L)      E10.4     Time for subschedule K of schedule L

         ACN(3,K,L)      E10.4     Value for subschedule K of schedule L
         K=1,KSUB(L)
         L=1,KCN

                               	    END OF INPUT DATA
                                            67

-------
                    APPENDIX B. MAIN DATA FILES FOR EXAMPLE 2
EXAMPLE: 2A (INFILTRATION)
  100100210
1000  30  1   10   9  12  0  0   2  0
  0.0000   1.0000  2.0000  3.0000  4.0000   5.0000  6.0000
  7.0000   8.0000
  0.0000   1.0000  2.0000  3.0000  4.0000   5.0000  6.0000
  7.0000   8.0000  9.0000 10.0000  11.0000
  2   1.0000
  0.0000   0.0005100.0000  0.1500   1.0300120.0000
  2.0000   0.8320  2.6900  1.5900
  1.0000   0.1000  0.0010  1.0000
  0   1
  4.0000   3.5000  -1.0000
 10.0000  5.0000  0.3500  0.0500   0.2000  5.0000  2.8000
  930
  000
  1   6
  2   5
  3   4
  4   3
  5   2
 100  10
 101   9
 102   8
 103   7
 45  1
 54  1
 63  1
  3   0.0000  -1.0000
    0.005 -0.1000
  9999.0000  -0.1000
  2   0.0000   0.0000
  9999.0000   0.0000
  2   0.0000   1.0000
  9999.0000   1.0000
  2   0.0000   2.0000
  9999.0000   2.0000
  2   0.0000   3.0000
  9999.0000   3.0000
  2   0.0000   4.0000
  9999.0000   4.0000
  2   0.0000   0.5000
  9999.0000   0.5000
  2   0.0000   1.5000
  9999.0000   1.5000
  2   0.0000   2.5000
  9999.0000   2.5000
  2   0.0000   3.5000
  9999.0000   3.5000
  000000
.00000000 .00000000 .00000000 .00000000 .00000000
  0
                                               68

-------
EXAMPLE: 2B (REDISTRIBUTION)..
1 1 0
1000 30 1
0.0000 1.
7.0000 8.
0.0000 1.
7.0000 8.
0 0.0000
0.0000 0.
2.0000 0.
1.0000 0.
10.0000 5
930
000
1 6
2 5
3 4
4 3
5 2
100 10
101 9
102 8
103 7
45 1
54 1
63 1
3 0.0000
1 0
10
0000
0000
0000
0000

0005
8320
1000
.0000














1
9
2.

2.
9.

25
2.
0.
0














-0.1000
2 1
12 0
0000

0000
0000

.0000
6900
0010
.3500














0.
0 0
0 2
3.0000

3.0000
10.0000

0.1500
1.5900
1.0000
0.0500














0001 -0.

0
4.0000

4.0000
11.0000

1.0300


0.2000














1000
                                         5.0000  6.0000

                                         5.0000  6.0000


                                         120.0000


                                         5.0000  2.8000
0.0002 -0.1000
 2   0.0000   0.0000
  9999.0000  0.0000
 2   0.0000   1.0000
  9999.0000  1.0000
 2   0.0000   2.0000
  9999.0000  2.0000
 2   0.0000   3.0000
  9999.0000  3.0000
 2   0.0000   4.0000
  9999.0000  4.0000
 2   0.0000   0.5000
  9999.0000  0.5000
 2   0.0000   1.5000
  9999.0000  1.5000
 2   0.0000   2.5000
  9999.0000  2.5000
 2   0.0000   3.5000
  9999.0000  3.5000
 220810
0.2000000 0.0400000 1.0000000 2.00000
493.0000 0.2400000 .00000000 .00000000.00000000.00000000.00000000
1000000.0 0.0020000 .00000000 .00000000 .00000000 .00000000 .00000000
.00009400 .00011320 0.7629000 877000.0
.00000010 .00000010 .00000100 827000.0
 0  177.89100.7443000
1
2
3
4
5
6
4 1
4 1
4 1
4 1
4 1
4 1
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
                                               69

-------
7410000
8410000
2.0000000 .0000000
  9999.00 .0000000
                                             70

-------
                  APPENDIX C. MAIN OUPUT FILES FOR EXAMPLE 2
              Cl. OUTPUT FOR EXAMPLE 2A  (STAGE 1 - INFILTRATION)
  TITLE: EXAMPLE: 2A (INFILTRATION)
  NUMBER OF NODES           =  108
  NUMBER OF ELEMENTS        =  88
  OIL PHASE INDEX, IOIL         =   1
  GAS PHASE INDEX, IAIR        =   0
  NUMBER OF ACTIVE PHASES    =   2
  NUMBER OF ROWS            =   9
  NUMBER OF COLUMNS         =  12
  IRAD                      =   0
  IRES                      =   0
  IRUN                      =   1
  ITRN                      =   0
  IDIM                       =   2
  NO. OF SCHEDULES           =   10
  IFEMC                     =   1
  KFEMC                     =   0
  ISOLVE                     =   2

—TEMPORAL DATA	
  STARTING TIME              =  0.0000
  TIME STEP AT T=0            =  0.0005
  MAX. SIMULATION TIME        = 100.0000
  MAX. TIME STEP              =  0.1500
  T. INCREMENT FACTOR        =  1.0300
  PRINT OUT T. INTERVAL        =120.0000

—FLUID PROPERTIES DATA-—
  DENSITY RATIO (O/W)          =  0.8320
  VISCOSITY RATIO (O/W)        =  2.0000
  BAD                       =  2.6900
  BOW                      =  1.5900

—CONTROL DATA	
  TIME WEIGHTING FACTOR      =  1.0000
  ABSOLUTE CONVERGENCE     =  0.1000
  RELATIVE ERROR            = 0.00100
  UP. STREAM FACTOR FOR FLOW =    1.00
  MAX. CUM. INTAKE (TINF)       =   1.000
  JINF                       =      2
  MAX. ITERATIONS            =     30
  MAX. TIME STEPS            =   1000
  TOL1
  TOL2
  TOL3
                   =  0.1000
                   =  0.0002
                   = 0.00020
ELEMENT NO  	NODES CONNECTED—
     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
 1
 2
 3
 4
 5
 6
 7
 8
10
11
12
13
14
15
16
17
10
11
12
13
14
15
16
17
19
20
21
22
23
24
25
26
11
12
13
14
15
16
17
18
20
21
22
23
24
25
26
27
2
3
4
5
6
7
8
9
11
12
13
14
15
16
17
18
                                            71

-------
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
19
20
21
22
23
24
25
26
28
29
30
31
32
33
34
35
37
38
39
40
41
42
43
44
46
47
48
49
50
51
52
53
55
56
57
58
59
60
61
62
64
65
66
67
68
69
70
71
73
74
75
76
77
78
79
80
82
83
84
85
86
87
88
89
91
92
93
94
95
96
97
98
28
29
30
31
32
33
34
35
37
38
39
40
41
42
43
44
46
47
48
49
50
51
52
53
55
56
57
58
59
60
61
62
64
65
66
67
68
69
70
71
73
74
75
76
77
78
79
80
82
83
84
85
86
87
88
89
91
92
93
94
95
96
97
98
100
101
102
103
104
105
106
107
29
30
31
32
33
34
35
36
38
39
40
41
42
43
44
45
47
48
49
50
51
52
53
54
56
57
58
59
60
61
62
63
65
66
67
68
69
70
71
72
74
75
76
77
78
79
80
81
83
84
85
86
87
88
89
90
92
93
94
95
96
97
98
99
101
102
103
104
105
106
107
108
20
21
22
23
24
25
26
27
29
30
31
32
33
34
35
36
38
39
40
41
42
43
44
45
47
48
49
50
51
52
53
54
56
57
58
59
60
61
62
63
65
66
67
68
69
70
71
72
74
75
76
77
78
79
80
81
83
84
85
86
87
88
89
90
92
93
94
95
96
97
98
99
LINEAR INTERPOLATION COEFF FOR WATER
A1     A2     B
   4.0000000   3.5000000  -1.0000000
	MATERIAL PROPERTIES	
                                               72

-------
-MATERIAL TYPE-
KSWX
KRWY
POROSITY
SM
SROW
ALPHA
N
NUMBER OF TYPE 1 BC NODES (W)
NUMBER OF TYPE 1 BC NODES (O)
NUMBER OF TYPE 1 BC NODES (G)
IXBC3 0

NUMBER OF TYPE 2 BC ELEMENTS (W)
NUMBER OF TYPE 2 BC ELEMENTS (O)
NUMBER OF TYPE 2 BC ELEMENTS (G)
NODE SCHEDULE(W)
1 6
2 5
3 4
4 3
5 2
100 10
101 9
102 8
103 7
NODE SCHEDULE(O)
45 1
54 1
63 1
BOUNDARY CONDITION SCHEDULES
0.0000 -1.0000 ..SCHEDULE 1
0.0050 -0.1000
9999.0000 -0.1000
0.0000 0.0000 ..SCHEDULE 2
9999.0000 0.0000
0.0000 1.0000 ..SCHEDULE 3
9999.0000 1 .0000
0.0000 2.0000 ..SCHEDULE 4
9999.0000 2.0000
0.0000 3.0000 ..SCHEDULE 5
9999.0000 3.0000
0.0000 4.0000 ..SCHEDULE 6
9999.0000 4.0000
0.0000 0.5000 ..SCHEDULE 7
9999.0000 0.5000
0.0000 1.5000 ..SCHEDULE 8
9999.0000 1 .5000
0.1000E+02
0.5000E+01
0.3500
0.0500
0.2000
5.0000
2.8000
9
3
0

0
0
0






















                                            73

-------
  0.0000   2.5000  ..SCHEDULE 9
 9999.0000    2.5000
  0.0000   3.5000  ..SCHEDULE 10
 9999.0000    3.5000
'** END OF INPUT DATA ********
NODE
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
PROP X
1 0.0000000
1 0.0000000
1 0.0000000
1 0.0000000
1 0.0000000
1 0.0000000
1 0.0000000
1 0.0000000
1 0.0000000
1 1 .0000000
1 1 .0000000
1 1 .0000000
1 1 .0000000
1 1 .0000000
1 1 .0000000
1 1 .0000000
1 1 .0000000
1 1 .0000000
1 2.0000000
1 2.0000000
1 2.0000000
1 2.0000000
1 2.0000000
1 2.0000000
1 2.0000000
1 2.0000000
1 2.0000000
1 3.0000000
1 3.0000000
1 3.0000000
1 3.0000000
1 3.0000000
1 3.0000000
1 3.0000000
1 3.0000000
1 3.0000000
1 4.0000000
1 4.0000000
1 4.0000000
1 4.0000000
1 4.0000000
1 4.0000000
1 4.0000000
1 4.0000000
1 4.0000000
1 5.0000000
1 5.0000000
1 5.0000000
1 5.0000000
1 5.0000000
1 5.0000000
1 5.0000000
1 5.0000000
1 5.0000000
1 6.0000000
1 6.0000000
1 6.0000000
1 6.0000000
1 6.0000000
1 6.0000000
1 6.0000000
1 6.0000000
1 6.0000000
1 7.0000000
1 7.0000000
1 7.0000000
z
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
                                                  74

-------
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
1 7.0000000
1 7.0000000
1 7.0000000
1 7.0000000
1 7.0000000
1 7.0000000
1 8.0000000
1 8.0000000
1 8.0000000
1 8.0000000
1 8.0000000
1 8.0000000
1 8.0000000
1 8.0000000
1 8.0000000
1 9.0000000
1 9.0000000
1 9.0000000
1 9.0000000
1 9.0000000
1 9.0000000
1 9.0000000
1 9.0000000
1 9.0000000
1 10.0000000
1 10.0000000
1 10.0000000
1 10.0000000
1 10.0000000
1 10.0000000
1 10.0000000
1 10.0000000
1 10.0000000
1 1 1 .0000000
1 1 1 .0000000
1 1 1 .0000000
1 1 1 .0000000
1 1 1 .0000000
1 1 1 .0000000
1 1 1 .0000000
1 1 1 .0000000
1 1 1 .0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
HALF BANDWIDTH=
                    22
 NODE	HEADS - - -| 	SATURATIONS -
    HW    HO    HA   SW   SO    SA
.... INITIAL HEADS & SATURATIONS	
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
4.0000
3.0000
2.0000
1 .0000
0.0000
-1 .0000
-2.0000
-3.0000
-4.0000
3.9545
2.9545
1 .9545
0.9545
-0.0455
-1 .0455
-2.0455
-3.0455
-4.0455
3.9091
2.9091
1 .9091
0.9091
-0.0909
-1 .0909
4.0000
3.0000
2.0000
1 .0000
0.0000
-0.3715
-0.7430
-1.1145
-1 .4860
3.9545
2.9545
1 .9545
0.9545
-0.0169
-0.3884
-0.7599
-1.1314
-1 .5029
3.9091
2.9091
1 .9091
0.9091
-0.0338
-0.4053
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
1 .0000
1 .0000
1 .0000
1 .0000
1 .0000
0.1021
0.0651
0.0573
0.0543
1 .0000
1 .0000
1 .0000
1 .0000
0.9905
0.0982
0.0645
0.0571
0.0542
1 .0000
1 .0000
1 .0000
1 .0000
0.9385
0.0946
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.8979
0.9349
0.9427
0.9457
0.0000
0.0000
0.0000
0.0000
0.0095
0.9018
0.9355
0.9429
0.9458
0.0000
0.0000
0.0000
0.0000
0.0615
0.9054
                                             75

-------
25  -2.0909 -0.7768  0.0000  0.0639   0.0000   0.9361
26  -3.0909 -1.1483  0.0000  0.0569   0.0000   0.9431
27  -4.0909 -1.5198  0.0000  0.0542   0.0000   0.9458
28   3.8636  3.8636  0.0000  1.0000   0.0000   0.0000
29   2.8636  2.8636  0.0000  1.0000   0.0000   0.0000
30   1.8636  1.8636  0.0000  1.0000   0.0000   0.0000
31   0.8636  0.8636  0.0000  1.0000   0.0000   0.0000
32  -0.1364 -0.0507  0.0000  0.8365   0.0000   0.1635
33  -1.1364 -0.4222  0.0000  0.0915   0.0000   0.9085
34  -2.1364 -0.7936  0.0000  0.0634   0.0000   0.9366
35  -3.1364 -1.1651  0.0000  0.0567   0.0000   0.9433
36  -4.1364 -1.5366  0.0000  0.0541   0.0000   0.9459
37   3.8182  3.8182  0.0000  1.0000   0.0000   0.0000
38   2.8182  2.8182  0.0000  1.0000   0.0000   0.0000
39   1.8182  1.8182  0.0000  1.0000   0.0000   0.0000
40   0.8182  0.8182  0.0000  1.0000   0.0000   0.0000
41  -0.1818 -0.0675  0.0000  0.7095   0.0000   0.2905
42  -1.1818 -0.4390  0.0000  0.0887   0.0000   0.9113
43  -2.1818 -0.8105  0.0000  0.0629   0.0000   0.9371
44  -3.1818 -1.1820  0.0000  0.0565   0.0000   0.9435
45  -4.1818 -1.5535  0.0000  0.0540   0.0000   0.9460
46   3.7727  3.7727  0.0000  1.0000   0.0000   0.0000
47   2.7727  2.7727  0.0000  1.0000   0.0000   0.0000
48   1.7727  1.7727  0.0000  1.0000   0.0000   0.0000
49   0.7727  0.7727  0.0000  1.0000   0.0000   0.0000
50  -0.2273 -0.0844  0.0000  0.5872   0.0000   0.4128
51  -1.2273 -0.4559  0.0000  0.0862   0.0000   0.9138
52  -2.2273 -0.8274  0.0000  0.0624   0.0000   0.9376
53  -3.2273 -1.1989  0.0000  0.0564   0.0000   0.9436
54  -4.2273 -1.5704  0.0000  0.0539   0.0000   0.9461
55   3.7273  3.7273  0.0000  1.0000   0.0000   0.0000
56   2.7273  2.7273  0.0000  1.0000   0.0000   0.0000
57   1.7273  1.7273  0.0000  1.0000   0.0000   0.0000
58   0.7273  0.7273  0.0000  1.0000   0.0000   0.0000
59  -0.2727 -0.1013  0.0000  0.4843   0.0000   0.5157
60  -1.2727 -0.4728  0.0000  0.0839   0.0000   0.9161
61  -2.2727 -0.8443  0.0000  0.0620   0.0000   0.9380
62  -3.2727 -1.2158  0.0000  0.0562   0.0000   0.9438
63  -4.2727 -1.5873  0.0000  0.0538   0.0000   0.9462
64   3.6818  3.6818  0.0000  1.0000   0.0000   0.0000
65   2.6818  2.6818  0.0000  1.0000   0.0000   0.0000
66   1.6818  1.6818  0.0000  1.0000   0.0000   0.0000
67   0.6818  0.6818  0.0000  1.0000   0.0000   0.0000
68  -0.3182 -0.1182  0.0000  0.4031   0.0000   0.5969
69  -1.3182 -0.4897  0.0000  0.0818   0.0000   0.9182
70  -2.3182 -0.8612  0.0000  0.0615   0.0000   0.9385
71  -3.3182 -1.2327  0.0000  0.0561   0.0000   0.9439
72  -4.3182 -1.6042  0.0000  0.0538   0.0000   0.9462
73   3.6364  3.6364  0.0000  1.0000   0.0000   0.0000
74   2.6364  2.6364  0.0000  1.0000   0.0000   0.0000
75   1.6364  1.6364  0.0000  1.0000   0.0000   0.0000
76   0.6364  0.6364  0.0000  1.0000   0.0000   0.0000
77  -0.3636 -0.1351  0.0000  0.3403   0.0000   0.6597
78  -1.3636 -0.5066  0.0000  0.0799   0.0000   0.9201
79  -2.3636 -0.8781  0.0000  0.0612   0.0000   0.9388
80  -3.3636 -1.2496  0.0000  0.0559   0.0000   0.9441
81  -4.3636 -1.6211  0.0000  0.0537   0.0000   0.9463
82   3.5909  3.5909  0.0000  1.0000   0.0000   0.0000
83   2.5909  2.5909  0.0000  1.0000   0.0000   0.0000
84   1.5909  1.5909  0.0000  1.0000   0.0000   0.0000
85   0.5909  0.5909  0.0000  1.0000   0.0000   0.0000
86  -0.4091 -0.1520  0.0000  0.2918   0.0000   0.7082
87  -1.4091 -0.5235  0.0000  0.0782   0.0000   0.9218
88  -2.4091 -0.8950  0.0000  0.0608   0.0000   0.9392
89  -3.4091 -1.2665  0.0000  0.0558   0.0000   0.9442
90  -4.4091 -1.6380  0.0000  0.0536   0.0000   0.9464
91   3.5455  3.5455  0.0000  1.0000   0.0000   0.0000
92   2.5455  2.5455  0.0000  1.0000   0.0000   0.0000
93   1.5455  1.5455  0.0000  1.0000   0.0000   0.0000
94   0.5455  0.5455  0.0000  1.0000   0.0000   0.0000
95  -0.4545 -0.1689  0.0000  0.2540   0.0000   0.7460
96  -1.4545 -0.5404  0.0000  0.0767   0.0000   0.9233
97  -2.4545 -0.9119  0.0000  0.0604   0.0000   0.9396
98  -3.4545 -1.2833  0.0000  0.0556   0.0000   0.9444
99  -4.4545 -1.6548  0.0000  0.0536   0.0000   0.9464
100  3.5000  3.5000  0.0000   1.0000   0.0000  0.0000
101   2.5000  2.5000  0.0000   1.0000   0.0000  0.0000
102  1.5000  1.5000  0.0000   1.0000   0.0000  0.0000
103  0.5000  0.5000  0.0000   1.0000   0.0000  0.0000
104  -0.5000  -0.1857   0.0000   0.2243  0.0000  0.7757
105  -1.5000  -0.5572   0.0000   0.0752  0.0000  0.9248
                                                        76

-------
106 -2.5000  -0.9287  0.0000  0.0601  0.0000  0.9399
107 -3.5000  -1.3002  0.0000  0.0555  0.0000  0.9445
108 -4.5000  -1.6717  0.0000  0.0535  0.0000  0.9465
  INITIAL VOLUME OF WATER =
  INITIAL VOLUME OF OIL =
  16.6242
0.0000
     ' END OF INITIAL CONDITIONS '
   FINAL OUTPUT AT TIME=  7.688650


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
QW-X
0.4221 E+00
0.4185E+00
0.4062E+00
0.3802E+00
0.3470E+00
-0.7999E-06
0.4426E-06
0.401 4E-06
0.3052E-06
0.4221 E+00
0.4185E+00
0.4062E+00
0.3802E+00
0.3470E+00
-0.7999E-06
0.4426E-06
0.401 4E-06
0.3052E-06
0.4253E+00
0.4224E+00
0.4126E+00
0.391 OE+00
0.3143E+00
-0.9549E-06
-0.1116E-06
-0.2646E-05
-0.4347E-05
0.431 3E+00
0.4294E+00
0.4237E+00
0.4127E+00
0.2534E+00
-0.1166E-05
-0.3291 E-05
-0.601 4E-05
-0.6509E-05
0.4390E+00
0.4383E+00
0.4370E+00
0.4386E+00
0.1798E+00
-0.1 297 E-05
-0.1 506 E-05
-0.5634E-06
0.1064E-05
0.4476E+00
0.4478E+00
0.4496E+00
0.4588E+00
0.1112E+00
-0.1 275 E-05
-0.1645E-06
0.1547E-06
0.2737E-06
0.4561 E+00
0.4568E+00
0.4599E+00
0.4689E+00
QW-Y
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1597E+01
-0.1584E-05
-0.3245E-08
0.6060E-07
-0.1785E-02
-0.1785E-02
-0.6135E-02
-0.1299E-01
-0.1481E-01
0.1141E+01
-0.1 189 E-05
-0.2385E-07
0.1250E-07
-0.3245E-02
-0.3245E-02
-0.1103E-01
-0.2378E-01
-0.2828E-01
0.651 7E+00
-0.1 070 E-05
-0.1 291 E-05
-0.8379E-06
-0.4166E-02
-0.4166E-02
-0.1388E-01
-0.2922E-01
-0.341 6E-01
0.2691 E+00
-0.2535E-05
-0.2653E-05
-0.1 086 E-05
-0.4508E-02
-0.4508E-02
-0.1455E-01
-0.2831 E-01
-0.3042E-01
0.5011 E-01
-0.3096E-05
-0.21 81 E-05
-0.2721 E-06
-0.4401 E-02
-0.4401 E-02
-0.1363E-01
-0.2354E-01
-0.2122E-01
-0.701 4E-06
-0.2988E-05
-0.2022E-05
-0.2126E-06
-0.4047E-02
-0.4047E-02
-0.1210E-01
-0.1883E-01
	 i
QO-X
0.1017E-06
0.9894E-07
0.9633E-07
0.8595E-07
0.7822E-07
-0.6475E-07
0.8164E-07
0.7232E-07
0.4740E-07
0.1017E-06
0.9894E-07
0.9633E-07
0.8595E-07
0.7822E-07
-0.6475E-07
0.8164E-07
0.7232E-07
0.4740E-07
0.1930E-06
0.1908E-06
0.1870E-06
0.1772E-06
0.641 OE-07
-0.6489E-07
-0.1 900 E-06
-0.1 425 E-05
-0.221 7E-05
0.2167E-06
0.2158E-06
0.2129E-06
0.2060E-06
0.6540E-07
-0.8122E-07
-0.1632E-04
-0.241 9E-02
-0.1109E-03
0.2227E-06
0.2224E-06
0.221 4E-06
0.2168E-06
0.681 5E-07
-0.1684E-06
-0.7888E-03
-0.5834E-02
-0.31 91 E-01
0.2254E-06
0.2255E-06
0.2256E-06
0.2292E-06
0.7171E-07
-0.2007E-06
-0.1 080 E-02
-0.1816E-02
O.OOOOE+00
0.2277E-06
0.2278E-06
0.2290E-06
0.2341 E-06

QO-Y
0.2706E-06
0.2706E-06
0.3962E-06
0.401 6E-06
0.3545E-06
-0.8472E-06
-0.1 445 E-05
-0.1 152 E-05
-0.1124E-05
0.2692E-06
0.2692E-06
0.3949E-06
0.3964E-06
0.3506E-06
-0.9187E-06
-0.1 372 E-05
-0.1 156 E-05
-0.1 137 E-05
0.2681 E-06
0.2681 E-06
0.3930E-06
0.391 5E-06
0.2941 E-06
-0.9832E-06
-0.1 435 E-05
-0.1774E-05
-0.1 533 E-05
0.2676E-06
0.2676E-06
0.391 5E-06
0.3881 E-06
0.2238E-06
-0.1056E-05
-0.2086E-04
-0.1 860 E-02
-0.5249E-04
0.2675E-06
0.2675E-06
0.391 OE-06
0.3858E-06
0.1495E-06
-0.1 175 E-05
-0.2344E-02
-0.3393E-01
-0.6187E-01
0.2675E-06
0.2675E-06
0.391 1 E-06
0.3876E-06
0.7071 E-07
-0.1 311 E-05
-0.1 273 E-01
-0.6047E-01
-0.6092E-01
0.2676E-06
0.2676E-06
0.391 6E-06
0.3902E-06

QA-X
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00

QA-Y
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
                                                    77

-------
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
NOC

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
0.5852E-01
-0.1146E-05
0.771 6E-06
0.4305E-06
0.3197E-06
0.4642E+00
0.4650E+00
0.4678E+00
0.471 9E+00
0.2650E-01
0.2622E-06
0.2115E-05
0.1164E-05
-0.4242E-06
0.471 3E+00
0.4723E+00
0.4743E+00
0.4745E+00
0.1044E-01
0.2448E-05
0.451 OE-05
0.6778E-05
0.7165E-05-
0.4773E+00
0.4789E+00
0.481 6E+00
0.4793E+00
0.3362E-02
0.3453E-05
0.1273E-05
0.431 1 E-05
0.6060E-05
0.481 8E+00
0.4838E+00
0.491 5E+00
0.4892E+00
-0.4673E-03
0.2195E-05
0.4639E-06
0.5026E-06
0.6111E-06
0.4841 E+00
0.4868E+00
0.4958E+00
0.5307E+00
-0.941 1 E-02
0.831 7E-06
0.4422E-06
0.4196E-06
0.3872E-06
-0.1320E-01
-0.2572E-05
-0.2486E-05
-0.2192E-05
-0.2681 E-06
-0.3625E-02
-0.3625E-02
-0.1073E-01
-0.1652E-01
-0.8860E-02
-0.3694E-05
-0.1 455 E-05
-0.2668E-05
-0.1 062 E-05
-0.3133E-02
-0.3133E-02
-0.971 5E-02
-0.1617E-01
-0.7672E-02
-0.2563E-05
0.2496E-05
-0.1534E-05
0.8685E-06
-0.2337E-02
-0.2337E-02
-0.8368E-02
-0.1701E-01
-0.891 5E-02
-0.7527E-06
0.2302E-05
-0.1503E-07
0.5821 E-08
-0.1336E-02
-0.1336E-02
-0.451 1 E-02
-0.1782E-01
-0.1411E-01
0.2502E-02
0.1256E-05
0.4344E-08
0.6004E-07
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
-0.3698E-01
0.1581E-01
0.1020E-05
-0.6946E-08
0.4384E-07
E i 	 HEADS 	 I I 	

HW
4.0000
3.0000
2.0000
1 .0000
0.0000
-1 .3201
-2.0033
-3.0026
-4.01 47
3.9578
2.9582
1 .9594
0.9620
-0.0348
-1 .2852
-2.0475
-3.0427
-4.0452
3.9153
2.9159
1.9181
0.9229
-0.0693
HO
3.8999
2.9597
1 .9692
0.9765
0.0027
-0.4904
-0.7443
-1.1157
-1 .4979
3.8796
2.9399
1 .9499
0.9593
-0.0129
-0.4774
-0.7607
-1.1301
-1 .5074
3.8410
2.9017
1.9125
0.9239
-0.0257
0.7283E-07
0.1106E-06
0.1172E-02
0.1 905 E-02
O.OOOOE+00
0.2297E-06
0.2299E-06
0.2311 E-06
0.2357E-06
0.6741 E-07
0.8794E-07
0.7176E-03
0.5653E-02
0.3206E-01
0.231 3E-06
0.231 6E-06
0.2329E-06
0.2367E-06
0.5157E-07
0.9286E-07
0.1763E-04
0.2490E-02
0.1123E-03
0.2308E-06
0.231 6E-06
0.2343E-06
0.2392E-06
0.2720E-07
0.1 480 E-06
0.4821 E-06
0.1 966 E-05
0.2781 E-05
0.2160E-06
0.2224E-06
0.2347E-06
0.2446E-06
-0.4665E-08
0.1442E-06
0.8648E-07
0.9557E-07
0.1274E-06
0.1000E-06
0.1926E-06
0.2306E-06
0.2632E-06
-0.5547E-07
0.8086E-07
0.8188E-07
0.7725E-07
0.7031 E-07
SATURATIONS
HA
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
-0.9907E-08
-0.1 292 E-05
-0.2035E-02
-0.3309E-01
-0.6193E-01
0.2677E-06
0.2677E-06
0.3922E-06
0.3925E-06
-0.9406E-07
-0.1 282 E-05
-0.1812E-04
-0.1 857 E-02
-0.5166E-04
0.2679E-06
0.2679E-06
0.3929E-06
0.3944E-06
-0.1 866 E-06
-0.1 261 E-05
-0.1 167 E-05
-0.1 897 E-05
-0.1 543 E-05
0.2683E-06
0.2683E-06
0.3942E-06
0.3969E-06
-0.2926E-06
-0.1 201 E-05
-0.1 000 E-05
-0.1 155 E-05
-0.1 136 E-05
0.271 5E-06
0.271 5E-06
0.4003E-06
0.401 8E-06
-0.4172E-06
-0.1 126 E-05
-0.1 029 E-05
-0.1 150 E-05
-0.1 120 E-05
0.3178E-06
0.3178E-06
0.4193E-06
0.41 81 E-06
-0.5765E-06
-0.1058E-05
-0.1029E-05
-0.1153E-05
-0.1123E-05


SW
1 .0000
1 .0000
1 .0000
1 .0000
1 .0000
0.0817
0.0650
0.0573
0.0543
1 .0000
1 .0000
1 .0000
1 .0000
0.9952
0.0833
0.0644
0.0571
0.0543
1 .0000
1 .0000
1 .0000
1 .0000
0.9696
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00

SO
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0003
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0004
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00

SA
0.0000
0.0000
0.0000
0.0000
0.0000
0.9183
0.9350
0.9427
0.9457
0.0000
0.0000
0.0000
0.0000
0.0045
0.9167
0.9356
0.9429
0.9457
0.0000
0.0000
0.0000
0.0000
0.0300
78

-------
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
-1 .2503
-2.0364
-2.7781
-3.6106
3.8721
2.8730
1 .8757
0.8816
-0.1045
-1.2143
-1 .7073
-2.1767
-2.9596
3.8282
2.8291
1 .8320
0.8377
-0.1412
-1.1760
-1 .5567
-2.1204
-3.0660
3.7835
2.7844
1.7871
0.7918
-0.1798
-1.1379
-1 .5402
-2.1359
-3.0933
3.7379
2.7387
1.7411
0.7449
-0.2190
-1.1146
-1.6174
-2.1789
-3.1253
3.6915
2.6922
1 .6943
0.6976
-0.2553
-1.1199
-1 .8289
-2.2953
-3.0829
3.6443
2.6449
1 .6469
0.6501
-0.2830
-1.1648
-2.2799
-2.9731
-3.7994
3.5966
2.5971
1 .5987
0.6021
-0.2977
-1 .2444
-2.4073
-3.4043
-4.4054
3.5484
2.5487
1 .5496
0.5531
-0.2952
-1 .3221
-2.4537
-3.4545
-4.4665
3.5000
2.5000
1 .5000
0.5000
-0.2653
-0.4645
-0.7227
-0.8451
-1 .0640
3.7976
2.8586
1 .8699
0.8827
-0.0388
-0.4482
-0.3644
-0.1919
-0.3182
3.7531
2.8141
1 .8257
0.8393
-0.0524
-0.4145
-0.1859
-0.1142
-0.1000
3.7080
2.7690
1 .7806
0.7935
-0.0668
-0.3744
-0.1407
-0.1011
-0.1000
3.6625
2.7234
1 .7348
0.7467
-0.0813
-0.3965
-0.1897
-0.1149
-0.1000
3.6165
2.6774
1 .6885
0.6995
-0.0948
-0.4141
-0.3729
-0.1923
-0.3192
3.5703
2.6311
1 .6420
0.6522
-0.1051
-0.4327
-0.7978
-0.871 1
-1 .0857
3.5241
2.5848
1.5951
0.6044
-0.1106
-0.4623
-0.8942
-1 .2642
-1.6419
3.4809
2.5403
1 .5482
0.5555
-0.1097
-0.491 1
-0.9115
-1 .2834
-1 .6674
3.4609
2.5018
1.5021
0.5028
-0.0986
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0850
0.0639
0.0569
0.0542
1 .0000
1 .0000
1 .0000
1 .0000
0.9121
0.0866
0.0634
0.0566
0.0540
1 .0000
1 .0000
1 .0000
1 .0000
0.8233
0.0870
0.0629
0.0565
0.0532
1 .0000
1 .0000
1 .0000
1 .0000
0.7149
0.0868
0.0624
0.0563
0.0532
1 .0000
1 .0000
1 .0000
1 .0000
0.6078
0.0911
0.0620
0.0562
0.0531
1 .0000
1 .0000
1 .0000
1 .0000
0.5208
0.0923
0.0616
0.0560
0.0537
1 .0000
1 .0000
1 .0000
1 .0000
0.4638
0.0897
0.0612
0.0560
0.0538
1 .0000
1 .0000
1 .0000
1 .0000
0.4369
0.0853
0.0608
0.0558
0.0537
1 .0000
1 .0000
1 .0000
1 .0000
0.4414
0.0816
0.0604
0.0556
0.0536
1 .0000
1 .0000
1 .0000
1 .0000
0.4994
0.0000
0.0019
0.0050
0.0037
0.0000
0.0000
0.0000
0.0000
0.0005
0.0007
0.0406
0.1584
0.0646
0.0000
0.0000
0.0000
0.0000
0.0004
0.0059
0.1611
0.3639
0.4383
0.0000
0.0000
0.0000
0.0000
0.0004
0.0146
0.2604
0.4290
0.4384
0.0000
0.0000
0.0000
0.0000
0.0005
0.0053
0.1562
0.3614
0.4384
0.0000
0.0000
0.0000
0.0000
0.0004
0.0006
0.0402
0.1584
0.0646
0.0000
0.0000
0.0000
0.0000
0.0003
0.0000
0.0021
0.0053
0.0038
0.0000
0.0000
0.0000
0.0000
0.0003
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0003
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0002
0.9150
0.9342
0.9381
0.9421
0.0000
0.0000
0.0000
0.0000
0.0874
0.9127
0.8961
0.7850
0.8814
0.0000
0.0000
0.0000
0.0000
0.1763
0.9071
0.7760
0.5796
0.5085
0.0000
0.0000
0.0000
0.0000
0.2847
0.8986
0.6772
0.5147
0.5085
0.0000
0.0000
0.0000
0.0000
0.3917
0.9036
0.7818
0.5825
0.5085
0.0000
0.0000
0.0000
0.0000
0.4788
0.9071
0.8982
0.7856
0.8818
0.0000
0.0000
0.0000
0.0000
0.5359
0.9103
0.9367
0.9387
0.9424
0.0000
0.0000
0.0000
0.0000
0.5628
0.9147
0.9392
0.9442
0.9463
0.0000
0.0000
0.0000
0.0000
0.5583
0.9183
0.9396
0.9444
0.9464
0.0000
0.0000
0.0000
0.0000
0.5004
79

-------
105      -1.3656        -0.5073      0.0000       0.0798       0.0001       0.9201
106      -2.4979        -0.9279      0.0000       0.0601       0.0000       0.9399
107      -3.4965        -1.2988      0.0000       0.0555       0.0000       0.9445
108      -4.5053        -1.6815      0.0000       0.0535       0.0000       0.9465

   FINAL VOLUME OF WATER =      17.0398
   FINAL VOLUME OF OIL  =       1.0153
TOTAL ITERATIONS     217
TOTAL TIME STEPS     211
                                                    80

-------
              OUTPUT FOR EXAMPLE 2B (STAGE 1 - REDISTRIBUTION )
  TITLE: EXAMPLE: 2B (REDISTRIBUTION)..
  NUMBER OF NODES            =  108
  NUMBER OF ELEMENTS         =   88
  OIL PHASE INDEX, IOIL          =    1
  GAS PHASE INDEX, IAIR         =    0
  NUMBER OF ACTIVE PHASES     =    2
  NUMBER OF ROWS            =    9
  NUMBER OF COLUMNS         =   12
  IRAD                       =    0
  IRES                       =    1
  IRUN                       =    1
  ITRN                       =    1
  IDIM                        =    2
  NO. OF SCHEDULES            =   10
  IFEMC                      =    1
  KFEMC                      =    0
  ISOLVE                      =    2

—TEMPORAL DATA	
  STARTING TIME               =   0.0000
  TIME STEP AT T=0             =   0.0005
  MAX. SIMULATION TIME         =  25.0000
  MAX. TIME STEP               =   0.1500
  T. INCREMENT FACTOR         =   1.0300
  PRINT OUT T. INTERVAL         =120.0000

—FLUID PROPERTIES DATA-—
  DENSITY RATIO (O/W)          =   0.8320
  VISCOSITY RATIO (O/W)         =   2.0000
  BAD                        =   2.6900
  BOW                       =   1.5900

—CONTROL DATA	
  TIME WEIGHTING FACTOR       =   1.0000
  ABSOLUTE CONVERGENCE      =   0.1000
  RELATIVE ERROR             =  0.00100
  UP. STREAM FACTOR FOR FLOW  =     1.00
  MAX. CUM. INTAKE (TINF)        =    0.000
  JINF                        =       0
  MAX. ITERATIONS             =      30
  MAX. TIME STEPS             =    1000
  TOL1
  TOL2
  TOL3
                   =   0.1000
                   =   0.0002
                   =  0.00020
ELEMENT NO  	NODES CONNECTED—
     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
 1
 2
 3
 4
 5
 6
 7
 8
10
11
12
13
14
15
16
17
19
10
11
12
13
14
15
16
17
19
20
21
22
23
24
25
26
28
11
12
13
14
15
16
17
18
20
21
22
23
24
25
26
27
29
2
3
4
5
6
7
8
9
11
12
13
14
15
16
17
18
20
                                             81

-------
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
20
21
22
23
24
25
26
28
29
30
31
32
33
34
35
37
38
39
40
41
42
43
44
46
47
48
49
50
51
52
53
55
56
57
58
59
60
61
62
64
65
66
67
68
69
70
71
73
74
75
76
77
78
79
80
82
83
84
85
86
87
88
89
91
92
93
94
95
96
97
98
29
30
31
32
33
34
35
37
38
39
40
41
42
43
44
46
47
48
49
50
51
52
53
55
56
57
58
59
60
61
62
64
65
66
67
68
69
70
71
73
74
75
76
77
78
79
80
82
83
84
85
86
87
88
89
91
92
93
94
95
96
97
98
100
101
102
103
104
105
106
107
30
31
32
33
34
35
36
38
39
40
41
42
43
44
45
47
48
49
50
51
52
53
54
56
57
58
59
60
61
62
63
65
66
67
68
69
70
71
72
74
75
76
77
78
79
80
81
83
84
85
86
87
88
89
90
92
93
94
95
96
97
98
99
101
102
103
104
105
106
107
108
21
22
23
24
25
26
27
29
30
31
32
33
34
35
36
38
39
40
41
42
43
44
45
47
48
49
50
51
52
53
54
56
57
58
59
60
61
62
63
65
66
67
68
69
70
71
72
74
75
76
77
78
79
80
81
83
84
85
86
87
88
89
90
92
93
94
95
96
97
98
99
	MATERIAL PROPERTIES	
-MATERIAL TYPE--
  KSWX
  KRWY
  POROSITY
= 0.1000E+02
= 0.5000E+01
= 0.3500
                                               82

-------
SM
SROW
ALPHA
N
NUMBER OF TYPE 1 BC NODES (W)
NUMBER OF TYPE 1 BC NODES (O)
NUMBER OF TYPE 1 BC NODES (G)
IXBC3 0

NUMBER OF TYPE 2 BC ELEMENTS (W)
NUMBER OF TYPE 2 BC ELEMENTS (O)
NUMBER OF TYPE 2 BC ELEMENTS (G)
NODE SCHEDULE(W)
1 6
2 5
3 4
4 3
5 2
100 10
101 9
102 8
103 7
NODE SCHEDULE(O)
45 1
54 1
63 1
BOUNDARY CONDITION SCHEDULES
0.0000 -0.1000 ..SCHEDULE 1
0.0001 -0.1000
0.0002 -0.1000
0.0000 0.0000 ..SCHEDULE 2
9999.0000 0.0000
0.0000 1.0000 ..SCHEDULE 3
9999.0000 1 .0000
0.0000 2.0000 ..SCHEDULE 4
9999.0000 2.0000
0.0000 3.0000 ..SCHEDULE 5
9999.0000 3.0000
0.0000 4.0000 ..SCHEDULE 6
9999.0000 4.0000
0.0000 0.5000 ..SCHEDULE 7
9999.0000 0.5000
0.0000 1.5000 ..SCHEDULE 8
9999.0000 1 .5000
0.0000 2.5000 ..SCHEDULE 9
9999.0000 2.5000
0.0000 3.5000 ..SCHEDULE 10
9999.0000 3.5000
0.0500
0.2000
5.0000
2.8000
9
3
0

0
0
0


























83

-------
—TRANSPORT DATA	
  KUNIT                    =  2
  NO. OF SPECIES                 =  2
  NO. OF NODES WITH TYPE 1 BC           =  0
  NO. OF ELEMENTS WITH TYPE 2 BC          =8
  NO. OF SCHEDULES DEFINING BC          =  1
  EVAPORATIVE BC INDEX          =         0
  LONGITUDINAL DISPERSIVITY       =      0.2000
  TRANSVERSE DISPERSIVITY        =      0.0400
  DENSITY OF PURE WATER          =  0.1000E+07
  DENSITY OF UNCONTAMINATED AIR  =  0.1200E+04
  UPSTREAM WEIGHT               =      1.0000
  TIME TO START UPDATE           =      2.0000

—SPECIES DATA-—


—SPECIES NO:	
  EQUILIBRIUM PARTITION COEFF. BET. OIL AND WATER     = 493.0000
  EQUILIBRIUM PARTITION COEFF. BET. AIR AND WATER     =   0.2400
  EQUILIBRIUM PARTITION COEFF. BET. SOILD AND WATER   =   0.0000
  FIRST ORDER DECAY COEFF. IN WATER                 =   0.0000
  FIRST ORDER DECAY COEFF. IN OIL                    =   0.0000
  FIRST ORDER DECAY COEFF. IN AIR                    =   0.0000
  FIRST ORDER DECAY COEFF. ON SOLID                 =   0.0000
  DIFFUSION COEFF. IN WATER                         =   0.0001
  DIFFUSION COEFF. IN OIL                             =   0.0001
  DIFFUSION COEFF. IN AIR                             =   0.7629
  LIQUID DENSITY OF PURE SPECIES                     =   **********

—SPECIES NO:	
  EQUILIBRIUM PARTITION COEFF. BET. OIL AND WATER      =   **********
  EQUILIBRIUM PARTITION COEFF. BET. AIR AND WATER     =    0.0020
  EQUILIBRIUM PARTITION COEFF. BET. SOILD AND WATER   =    0.0000
  FIRST ORDER DECAY COEFF. IN WATER                 =    0.0000
  FIRST ORDER DECAY COEFF. IN OIL                    =    0.0000
  FIRST ORDER DECAY COEFF. IN AIR                    =    0.0000
  FIRST ORDER DECAY COEFF. ON SOLID                 =    0.0000
  DIFFUSION COEFF. IN WATER                         =    0.0000
  DIFFUSION COEFF. IN OIL                             =    0.0000
  DIFFUSION COEFF. IN AIR                             =    0.0000
  LIQUID DENSITY OF PURE SPECIES                     =   **********
 ELEMENT|-SIDE |	SCHEDULE FOR SPECIES NO:-
        12345

  141     0
  241     0
  341     0
  441     0
  541     0
  641     0
  741     0
  841     0

TRANSPORT BOUNDARY CONDITION SCHEDULES
   0.0000   0.0000  ..SCHEDULE 1

  9999.0000    0.0000

**** END OF INPUT DATA ********
                                               84

-------
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
1 0.0000000
1 0.0000000
1 0.0000000
1 0.0000000
1 0.0000000
1 0.0000000
1 0.0000000
1 0.0000000
1 0.0000000
1 1 .0000000
1 1 .0000000
1 1 .0000000
1 1 .0000000
1 1 .0000000
1 1 .0000000
1 1 .0000000
1 1 .0000000
1 1 .0000000
1 2.0000000
1 2.0000000
1 2.0000000
1 2.0000000
1 2.0000000
1 2.0000000
1 2.0000000
1 2.0000000
1 2.0000000
1 3.0000000
1 3.0000000
1 3.0000000
1 3.0000000
1 3.0000000
1 3.0000000
1 3.0000000
1 3.0000000
1 3.0000000
1 4.0000000
1 4.0000000
1 4.0000000
1 4.0000000
1 4.0000000
1 4.0000000
1 4.0000000
1 4.0000000
1 4.0000000
1 5.0000000
1 5.0000000
1 5.0000000
1 5.0000000
1 5.0000000
1 5.0000000
1 5.0000000
1 5.0000000
1 5.0000000
1 6.0000000
1 6.0000000
1 6.0000000
1 6.0000000
1 6.0000000
1 6.0000000
1 6.0000000
1 6.0000000
1 6.0000000
1 7.0000000
1 7.0000000
1 7.0000000
1 7.0000000
1 7.0000000
1 7.0000000
1 7.0000000
1 7.0000000
1 7.0000000
1 8.0000000
1 8.0000000
1 8.0000000
1 8.0000000
1 8.0000000
1 8.0000000
1 8.0000000
1 8.0000000
1 8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
85

-------
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
1 9.0000000
1 9.0000000
1 9.0000000
1 9.0000000
1 9.0000000
1 9.0000000
1 9.0000000
1 9.0000000
1 9.0000000
1 10.0000000
1 10.0000000
1 10.0000000
1 10.0000000
1 10.0000000
1 10.0000000
1 10.0000000
1 10.0000000
1 10.0000000
1 1 1 .0000000
1 1 1 .0000000
1 1 1 .0000000
1 1 1 .0000000
1 1 1 .0000000
1 1 1 .0000000
1 1 1 .0000000
1 1 1 .0000000
1 1 1 .0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
0.0000000
1 .0000000
2.0000000
3.0000000
4.0000000
5.0000000
6.0000000
7.0000000
8.0000000
HALF BANDWIDTH=
                    22
 NODE	HEADS---|
     HW    HO     HA
--SATURATIONS	
 SW    SO    SA
.... INITIAL HEADS & SATURATIONS	
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
4.0000
3.0000
2.0000
1 .0000
0.0000
-1 .3200
-2.0030
-3.0030
-4.0150
3.9580
2.9580
1 .9590
0.9620
-0.0348
-1 .2850
-2.0480
-3.0430
-4.0450
3.9150
2.9160
1.9180
0.9229
-0.0693
-1 .2500
-2.0360
-2.7780
-3.6110
3.8720
2.8730
1 .8760
0.8816
-0.1045
-1.2140
-1 .7070
-2.1770
-2.9600
3.8280
2.8290
1 .8320
3.9000
2.9600
1 .9690
0.9765
0.0027
-0.4904
-0.7443
-1.1160
-1 .4980
3.8800
2.9400
1 .9500
0.9593
-0.0129
-0.4774
-0.7607
-1.1300
-1 .5070
3.8410
2.9020
1.9130
0.9239
-0.0257
-0.4645
-0.7227
-0.8451
-1 .0640
3.7980
2.8590
1 .8700
0.8827
-0.0388
-0.4482
-0.3644
-0.1919
-0.3182
3.7530
2.8140
1 .8260
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
1 .0000
1 .0000
1 .0000
1 .0000
1 .0000
0.0817
0.0650
0.0573
0.0543
1 .0000
1 .0000
1 .0000
1 .0000
0.9952
0.0833
0.0644
0.0571
0.0543
1 .0000
1 .0000
1 .0000
1 .0000
0.9696
0.0850
0.0639
0.0569
0.0542
1 .0000
1 .0000
1 .0000
1 .0000
0.9120
0.0866
0.0634
0.0566
0.0540
1 .0000
1 .0000
1 .0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0003
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0004
0.0000
0.0019
0.0050
0.0037
0.0000
0.0000
0.0000
0.0000
0.0006
0.0007
0.0406
0.1584
0.0647
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.9183
0.9350
0.9427
0.9457
0.0000
0.0000
0.0000
0.0000
0.0045
0.9167
0.9356
0.9429
0.9457
0.0000
0.0000
0.0000
0.0000
0.0300
0.9150
0.9342
0.9381
0.9421
0.0000
0.0000
0.0000
0.0000
0.0874
0.9127
0.8961
0.7850
0.8814
0.0000
0.0000
0.0000
                                             86

-------
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
0.8377
-0.1412
-1.1760
-1 .5570
-2.1200
-3.0660
3.7830
2.7840
1 .7870
0.7918
-0.1798
-1.1380
-1 .5400
-2.1360
-3.0930
3.7380
2.7390
1.7410
0.7449
-0.2190
-1.1150
-1.6170
-2.1790
-3.1250
3.6910
2.6920
1 .6940
0.6976
-0.2553
-1.1200
-1 .8290
-2.2950
-3.0830
3.6440
2.6450
1 .6470
0.6501
-0.2830
-1.1650
-2.2800
-2.9730
-3.7990
3.5970
2.5970
1 .5990
0.6021
-0.2977
-1 .2440
-2.4070
-3.4040
-4.4050
3.5480
2.5490
1 .5500
0.5531
-0.2952
-1 .3220
-2.4540
-3.4550
-4.4670
3.5000
2.5000
1 .5000
0.5000
-0.2653
-1 .3660
-2.4980
-3.4970
-4.5050
0.8393
-0.0524
-0.4145
-0.1859
-0.1142
-0.1000
3.7080
2.7690
1.7810
0.7935
-0.0668
-0.3744
-0.1407
-0.1011
-0.1000
3.6620
2.7230
1 .7350
0.7467
-0.0813
-0.3965
-0.1897
-0.1149
-0.1000
3.6170
2.6770
1 .6890
0.6995
-0.0948
-0.4141
-0.3729
-0.1923
-0.3192
3.5700
2.6310
1 .6420
0.6522
-0.1051
-0.4327
-0.7978
-0.871 1
-1 .0860
3.5240
2.5850
1 .5950
0.6044
-0.1106
-0.4623
-0.8942
-1 .2640
-1 .6420
3.4810
2.5400
1 .5480
0.5555
-0.1097
-0.491 1
-0.9115
-1 .2830
-1 .6670
3.4610
2.5020
1 .5020
0.5028
-0.0986
-0.5073
-0.9279
-1 .2990
-1.6810
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
1 .0000
0.8231
0.0870
0.0629
0.0565
0.0532
1 .0000
1 .0000
1 .0000
1 .0000
0.7147
0.0868
0.0624
0.0563
0.0532
1 .0000
1 .0000
1 .0000
1 .0000
0.6077
0.0911
0.0620
0.0562
0.0531
1 .0000
1 .0000
1 .0000
1 .0000
0.5207
0.0923
0.0616
0.0560
0.0537
1 .0000
1 .0000
1 .0000
1 .0000
0.4638
0.0897
0.0612
0.0560
0.0538
1 .0000
1 .0000
1 .0000
1 .0000
0.4369
0.0853
0.0608
0.0558
0.0537
1 .0000
1 .0000
1 .0000
1 .0000
0.4414
0.0816
0.0604
0.0556
0.0536
1 .0000
1 .0000
1 .0000
1 .0000
0.4994
0.0798
0.0601
0.0555
0.0535
INITIAL VOLUME OF WATER =
INITIAL VOLUME OF
OIL =
1.
0.0000
0.0006
0.0059
0.1612
0.3640
0.4383
0.0000
0.0000
0.0000
0.0000
0.0006
0.0146
0.2604
0.4292
0.4384
0.0000
0.0000
0.0000
0.0000
0.0005
0.0054
0.1563
0.3612
0.4384
0.0000
0.0000
0.0000
0.0000
0.0005
0.0006
0.0402
0.1585
0.0646
0.0000
0.0000
0.0000
0.0000
0.0005
0.0000
0.0021
0.0053
0.0038
0.0000
0.0000
0.0000
0.0000
0.0003
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0003
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0002
0.0001
0.0000
0.0000
0.0000
1 7.0396
0158
0.0000
0.1763
0.9071
0.7760
0.5795
0.5085
0.0000
0.0000
0.0000
0.0000
0.2847
0.8986
0.6772
0.5145
0.5085
0.0000
0.0000
0.0000
0.0000
0.3918
0.9036
0.7818
0.5826
0.5085
0.0000
0.0000
0.0000
0.0000
0.4788
0.9070
0.8982
0.7856
0.8818
0.0000
0.0000
0.0000
0.0000
0.5357
0.9103
0.9367
0.9387
0.9424
0.0000
0.0000
0.0000
0.0000
0.5628
0.9147
0.9392
0.9442
0.9463
0.0000
0.0000
0.0000
0.0000
0.5583
0.9183
0.9396
0.9444
0.9464
0.0000
0.0000
0.0000
0.0000
0.5004
0.9201
0.9399
0.9445
0.9465


—INITIAL CONC OF SPECIES
  NODE-I-CONCENTRATIONS IN PHASE-
     WATER  OIL   GAS   SOLID
                                             87

-------
SPECIES NUMBER =
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1779E+03
0.1779E+03
0.1779E+03
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1779E+03
0.1779E+03
0.1779E+03
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1779E+03
0.1779E+03
0.1779E+03
0.1779E+03
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1779E+03
0.1779E+03
0.1779E+03
0.1779E+03
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1779E+03
0.1779E+03
0.1779E+03
0.1779E+03
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1779E+03
0.1779E+03
0.1779E+03
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.8770E+05
0.8770E+05
0.8770E+05
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.8770E+05
0.8770E+05
0.8770E+05
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.8770E+05
0.8770E+05
0.8770E+05
0.8770E+05
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.8770E+05
0.8770E+05
0.8770E+05
0.8770E+05
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.8770E+05
0.8770E+05
0.8770E+05
0.8770E+05
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.8770E+05
0.8770E+05
0.8770E+05
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.4269E+02
0.4269E+02
0.4269E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.4269E+02
0.4269E+02
0.4269E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.4269E+02
0.4269E+02
0.4269E+02
0.4269E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.4269E+02
0.4269E+02
0.4269E+02
0.4269E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.4269E+02
0.4269E+02
0.4269E+02
0.4269E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.4269E+02
0.4269E+02
0.4269E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00

-------
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
O.OOOOE+00
0.1779E+03
0.1779E+03
0.1779E+03
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.8770E+05
0.8770E+05
0.8770E+05
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.4269E+02
0.4269E+02
0.4269E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
TOTAL MASS IN
  WATER    OIL     GAS     SOLID
0.7996E+02  0.8889E+05  0.2439E+03  O.OOOOE+00
SPECIES NUMBER = 2
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7443E+00
0.7443E+00
0.7443E+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7443E+00
0.7443E+00
0.7443E+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7443E+06
0.7443E+06
0.7443E+06
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7443E+06
0.7443E+06
0.7443E+06
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1489E-02
0.1489E-02
0.1489E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1489E-02
0.1489E-02
0.1489E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
                                                  89

-------
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7443E+00
0.7443E+00
0.7443E+00
0.7443E+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7443E+00
0.7443E+00
0.7443E+00
0.7443E+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7443E+00
0.7443E+00
0.7443E+00
0.7443E+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7443E+00
0.7443E+00
0.7443E+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7443E+00
0.7443E+00
0.7443E+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7443E+06
0.7443E+06
0.7443E+06
0.7443E+06
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7443E+06
0.7443E+06
0.7443E+06
0.7443E+06
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7443E+06
0.7443E+06
0.7443E+06
0.7443E+06
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7443E+06
0.7443E+06
0.7443E+06
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7443E+06
0.7443E+06
0.7443E+06
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1489E-02
0.1489E-02
0.1489E-02
0.1489E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1489E-02
0.1489E-02
0.1489E-02
0.1489E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1489E-02
0.1489E-02
0.1489E-02
0.1489E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1489E-02
0.1489E-02
0.1489E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1489E-02
0.1489E-02
0.1489E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
TOTAL MASS IN
  WATER   OIL     GAS     SOLID
0.3346E+00  0.7544E+06  0.8503E-02 O.OOOOE+00
                                                  90

-------
  ' END OF INITIAL CONDITIONS '
FINAL OUTPUT AT TIME= 25.088608


1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
QW-X QW-Y
0.4192E+00
0.4152E+00
0.4020E+00
0.3751 E+00
0.3409E+00
-0.9773E-06
0.3440E-06
-0.1079E-06
-0.1035E-05
0.4192E+00
0.4152E+00
0.4020E+00
0.3751 E+00
0.3409E+00
-0.9773E-06
0.3440E-06
-0.1079E-06
-0.1035E-05
0.4226E+00
0.4194E+00
0.4087E+00
0.3853E+00
0.3079E+00
-0.1046E-05
-0.2390E-05
-0.5598E-05
-0.7082E-05
0.4290E+00
0.4269E+00
0.4205E+00
0.4071 E+00
0.2449E+00
-0.4277E-06
-0.2284E-05
-0.2068E-05
-0.1743E-05
0.4374E+00
0.4365E+00
0.4346E+00
0.4359E+00
0.1 671 E+00
-0.3320E-06
-0.3724E-06
-0.2959E-06
0.1348E-05
0.4468E+00
0.4470E+00
0.4480E+00
0.4590E+00
0.971 7E-01
-0.341 4E-06
0.1697E-06
0.2380E-06
0.3180E-06
0.4563E+00
0.4571 E+00
0.4604E+00
0.4642E+00
0.4876E-01
-0.4432E-06
0.4063E-06
0.3625E-06
0.2823E-06
0.4651 E+00
0.4661 E+00
0.4704E+00
0.4668E+00
	 1



QO-X QO-Y QA-X QA-Y
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1367E+01
-0.1309E-05
-0.3773E-07
0.1448E-06
-0.2003E-02
-0.2003E-02
-0.6579E-02
-0.1348E-01
-0.1482E-01
0.941 7E+00
-0.9654E-06
-0.2637E-06
-0.3187E-06
-0.3604E-02
-0.3604E-02
-0.1192E-01
-0.251 8E-01
-0.291 4E-01
0.5021 E+00
-0.1992E-05
-0.1868E-05
-0.1060E-05
-0.4643E-02
-0.4643E-02
-0.1511E-01
-0.3179E-01
-0.3684E-01
0.1 761 E+00
-0.2966E-05
-0.1759E-05
-0.8981 E-06
-0.5096E-02
-0.5096E-02
-0.1603E-01
-0.3111E-01
-0.3387E-01
0.5902E-02
-0.2986E-05
-0.1721E-05
-0.7608E-07
-0.501 3E-02
-0.501 3E-02
-0.1550E-01
-0.2561 E-01
-0.2225E-01
-0.3263E-06
-0.2731 E-05
-0.1687E-05
-0.3607E-07
-0.4609E-02
-0.4609E-02
-0.1386E-01
-0.2293E-01
-0.1607E-01
-0.8627E-06
-0.2348E-05
-0.1 709 E-05
-0.761 8E-07
-0.4126E-02
-0.4126E-02
-0.1 171 E-01
-0.1737E-01
0.4340E-07
0.431 6E-07
0.4301 E-07
0.3029E-07
0.8596E-07
-0.6363E-07
0.5680E-07
-0.4951 E-07
-0.3251 E-06
0.4340E-07
0.431 6E-07
0.4301 E-07
0.3029E-07
0.8596E-07
-0.6363E-07
0.5680E-07
-0.4951 E-07
-0.3251 E-06
0.1194E-06
0.1205E-06
0.1224E-06
0.1024E-06
0.6280E-07
-0.1167E-06
-0.1 336 E-05
-0.4549E-05
-0.4242E-05
0.1726E-06
0.1717E-06
0.1 971 E-06
0.2066E-06
0.5863E-07
-0.2329E-04
-0.2776E-03
-0.4346E-03
-0.3233E-04
0.2024E-06
0.2046E-06
0.1977E-06
0.4797E-06
0.8589E-06
-0.7118E-03
-0.5119E-03-
-0.2024E-03
-0.1812E-03
0.2178E-06
0.2177E-06
0.2278E-06
0.8837E-07
0.681 OE-04
-0.4908E-03
-0.1255E-03
-0.3583E-04
-0.1 550 E-05
0.2244E-06
0.2224E-06
0.2466E-06
0.5066E-07
0.1052E-02
0.6265E-03
0.1505E-03
0.3932E-04
0.1 678 E-05
0.2250E-06
0.2274E-06
0.1478E-06
0.241 4E-05
0.1182E-06
0.1182E-06
0.2794E-06
0.3320E-06
0.2881 E-06
-0.8850E-06
-0.1 393 E-05
-0.1 158 E-05
-0.1064E-05
0.1 181 E-06
0.1 181 E-06
0.2793E-06
0.3257E-06
0.3159E-06
-0.9598E-06
-0.1 333 E-05
-0.1 211 E-05
-0.1 202 E-05
0.1186E-06
0.1186E-06
0.2803E-06
0.3157E-06
0.2961 E-06
-0.1 050 E-05
-0.1 942 E-05
-0.3073E-05
-0.1 775 E-05
0.1182E-06
0.1182E-06
0.2930E-06
0.3204E-06
0.2221 E-06
-0.4536E-04
-0.4721 E-03
-0.6820E-03
-0.4170E-04
0.1193E-06
0.1193E-06
0.2895E-06
0.461 4E-06
0.3427E-08
-0.2176E-02
0.3389E-02
-0.2178E-02
-0.6449E-03
0.1192E-06
0.1192E-06
0.2946E-06
0.391 7E-06
-0.2277E-03
-0.6604E-02
-0.491 4E-02
-0.2623E-02
-0.6587E-03
0.1182E-06
0.1182E-06
0.3067E-06
0.2937E-06
0.3446E-07
-0.1641E-02
-0.3180E-02
-0.2147E-02
-0.6443E-03
0.1194E-06
0.1194E-06
0.2668E-06
0.3567E-06
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
                                             91

-------
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
NO

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
0.2539E-01
-0.1014E-05
0.9755E-06
0.9241 E-06
-0.6285E-06
0.4730E+00
0.4742E+00
0.4766E+00
0.4726E+00
0.1118E-01
-0.1166E-05
0.2935E-05
0.2684E-05
0.2390E-05
0.4795E+00
0.481 3E+00
0.4838E+00
0.4791 E+00
0.4054E-02
0.2665E-06
0.3868E-05
0.7104E-05
0.8306E-05
0.4844E+00
0.4865E+00
0.4950E+00
0.4901 E+00
-0.2450E-03
0.1467E-05
0.6493E-06
0.1194E-05
0.2278E-05
0.4867E+00
0.4898E+00
0.4999E+00
0.5387E+00
-0.1119E-01
0.2720E-05
0.3748E-06
0.3878E-06
0.4176E-06
-0.1277E-01
-0.3569E-05
-0.1717E-05
-0.1735E-05
-0.8524E-06
-0.3526E-02
-0.3526E-02
-0.1050E-01
-0.1728E-01
-0.1149E-01
-0.1072E-04
-0.1 839 E-06
-0.1860E-05
-0.9990E-06
-0.2636E-02
-0.2636E-02
-0.9238E-02
-0.1889E-01
-0.1306E-01
-0.1202E-04
0.1656E-04
-0.2423E-06
-0.3982E-06
-0.1546E-02
-0.1546E-02
-0.5005E-02
-0.2045E-01
-0.1925E-01
-0.1018E-04
0.1716E-04
0.2996E-07
0.1440E-06
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
-0.481 1E-01
-0.6136E-05
0.1410E-04
0.3646E-07
0.1589E-06
0.5173E-04 -0
0.4759E-03 -0
0.4894E-03 -0
0.1999E-03 -0
0.1840E-03 -0
0.21 91 E-06 0
0.2186E-06 0
0.2583E-06 0
0.2527E-06 0
0.9638E-07 -0
0.1051E-04 -0
0.2705E-03 -0
0.4337E-03 -0
0.3142E-04 -0
0.1 951 E-06 0
0.2040E-06 0
0.2244E-06 0
0.2391 E-06 0
0.2621 E-07 -0
0.3494E-07 -0
0.1852E-05 -0
0.5422E-05 -0
0.4805E-05 -0
0.1342E-06 0
0.1556E-06 0
0.2054E-06 0
0.2448E-06 0
-0.1807E-08 -0
0.2836E-07 -0
0.1268E-06 -0
0.2593E-06 -0
0.5994E-06 -0
0.5587E-07 0
0.531 7E-08 0
0.1966E-06 0
0.2668E-06 0
-0.5346E-07 -0
0.5821 E-07 -0
1787 E-06
3702 E-04
4319E-03
6742 E-03
4026 E-04
11 91 E-06
11 91 E-06
2867E-06
3539E-06
1933 E-06
1389E-05
1742E-05
3133E-05
1777E-05
1236E-06
1236E-06
2969E-06
361 3E-06
2998 E-06
1385E-05
8335 E-06
1213E-05
1216E-05
1343E-06
1343E-06
321 8E-06
3809E-06
4231 E-06
1370E-05
7843 E-06
1146E-05
1046E-05
1 090E-06
1 090E-06
4175E-06
4160E-06
5832 E-06
1314E-05
0.6933E-07 -0.7787E-06
0.7376E-07 -0
0.8426E-07 -0
1 1 44E-05
1041E-05
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
DE 	 HEADS 	 I 	 SATURATIONS 	

HW HO
4.0000 3.7397
3.0000 2.8604
2.0000 1.9166
1.0000 0.9517
0.0000 0.0045
-1.2748 -0.4735
-2.0130 -0.7484
-3.0055 -1.1171
-4.0344 -1.5234
3.9581 3.7310
2.9585 2.8517
1 .9598 1 .9080
0.9625 0.9457
-0.0343 -0.0127
-1.2405 -0.4608
-2.0474 -0.7597
-2.9947 -1.1072
-3.9309 -1.4584
3.9158 3.7071
2.9165 2.8276
1.9189 1.8835
0.9240 0.9252
-0.0680 -0.0252
-1.2069 -0.4374
-1.8084 -0.4925
-2.4349 -0.4878
-3.2228 -0.6100
3.8729 3.6726
2.8739 2.7933
1 .8769 1 .8441
0.8832 0.8839
-0.1018 -0.0370
HA SW SO
0.0000 1 .0000
0.0000 1 .0000
0.0000 1 .0000
0.0000 1 .0000
0.0000 0.9999
0.0000 0.0837
0.0000 0.0649
0.0000 0.0572
0.0000 0.0543
0.0000 1 .0000
0.0000 1 .0000
0.0000 1 .0000
0.0000 1 .0000
0.0000 0.9954
0.0000 0.0853
0.0000 0.0644
0.0000 0.0572
0.0000 0.0545
0.0000 1 .0000
0.0000 1 .0000
0.0000 1 .0000
0.0000 1 .0000
0.0000 0.9701
0.0000 0.0862
0.0000 0.0639
0.0000 0.0569
0.0000 0.0540
0.0000 1 .0000
0.0000 1 .0000
0.0000 1 .0000
0.0000 1 .0000
0.0000 0.9133
SA
0.0000 0.0000
0.0000 0.0000
0.0000 0.0000
0.0000 0.0000
0.0001 0.0000
0.0001 0.9162
0.0000 0.9351
0.0000 0.9428
0.0000 0.9457
0.0000 0.0000
0.0000 0.0000
0.0000 0.0000
0.0000 0.0000
0.0003 0.0043
0.0001 0.9145
0.0000 0.9355
0.0001 0.9426
0.0000 0.9455
0.0000 0.0000
0.0000 0.0000
0.0000 0.0000
0.0000 0.0000
0.0013 0.0285
0.0028 0.9111
0.0176 0.9185
0.0252 0.9180
0.0174 0.9285
0.0000 0.0000
0.0000 0.0000
0.0000 0.0000
0.0000 0.0000
0.0092 0.0775



































































































92

-------
33  -1.1732 -0.2884  0.0000  0.0783  0.0534   0.8684
34  -1.5799 -0.2311  0.0000  0.0633  0.1069   0.8298
35  -2.2281  -0.2188  0.0000  0.0565  0.1255   0.8180
36  -3.0484 -0.3309  0.0000  0.0538  0.0603   0.8860
37   3.8292  3.6321  0.0000  1.0000  0.0000   0.0000
38   2.8302  2.7524  0.0000  1.0000  0.0000   0.0000
39   1.8334  1.8045  0.0000  1.0000  0.0000   0.0000
40   0.8396  0.7880  0.0000  1.0000  0.0000   0.0000
41  -0.1357 -0.0454  0.0000  0.8152  0.0571   0.1277
42  -1.1400 -0.1732  0.0000  0.0741  0.1716   0.7542
43  -1.5427 -0.1678  0.0000  0.0628  0.1933   0.7439
44  -2.1985 -0.1808  0.0000  0.0564  0.1759   0.7677
45  -3.1833 -0.2194  0.0000  0.0532  0.1282   0.8186
46   3.7845  3.5885  0.0000  1.0000  0.0000   0.0000
47   2.7855  2.7088  0.0000  1.0000  0.0000   0.0000
48   1.7886  1.7590  0.0000  1.0000  0.0000   0.0000
49   0.7937  0.7703  0.0000  1.0000  0.0000   0.0000
50  -0.1711  -0.0528  0.0000  0.6892  0.1318   0.1790
51  -1.1058 -0.1457  0.0000  0.0744  0.2339   0.6916
52  -1.5597 -0.1573  0.0000  0.0624  0.2164   0.7212
53  -2.2223 -0.1752  0.0000  0.0563  0.1858   0.7580
54  -3.2151  -0.2185  0.0000  0.0532  0.1292   0.8177
55   3.7389  3.5437  0.0000  1.0000  0.0000   0.0000
56   2.7398  2.6644  0.0000  1.0000  0.0000   0.0000
57   1.7426  1.7097  0.0000  1.0000  0.0000   0.0000
58   0.7473  0.7602  0.0000  0.9990  0.0010   0.0000
59  -0.2096 -0.0734  0.0000  0.6110  0.0542   0.3348
60  -1.0700 -0.1808  0.0000  0.0780  0.1543   0.7677
61  -1.6003 -0.1699  0.0000  0.0619  0.1901   0.7480
62  -2.2585 -0.1813  0.0000  0.0561  0.1753   0.7686
63  -3.2433 -0.2195  0.0000  0.0531  0.1282   0.8187
64   3.6924  3.4987  0.0000  1.0000  0.0000   0.0000
65   2.6932  2.6189  0.0000  1.0000  0.0000   0.0000
66   1.6955  1.6801  0.0000  1.0000  0.0000   0.0000
67   0.6990  0.7054  0.0000  0.9999  0.0001   0.0000
68  -0.2436 -0.0899  0.0000  0.5419  0.0093   0.4488
69  -1.0412 -0.2858  0.0000  0.0875  0.0454   0.8671
70  -1.6978 -0.2348  0.0000  0.0615  0.1055   0.8331
71  -2.3509 -0.2195  0.0000  0.0558  0.1255   0.8187
72  -3.1804 -0.3329  0.0000  0.0535  0.0599   0.8866
73   3.6451   3.4548  0.0000  1.0000  0.0000   0.0000
74   2.6458  2.5752  0.0000  1.0000  0.0000   0.0000
75   1.6479  1.6284  0.0000  1.0000  0.0000   0.0000
76   0.6513  0.6549  0.0000  1.0000  0.0000   0.0000
77  -0.2687 -0.0998  0.0000  0.4891  0.0034   0.5075
78  -1.0281  -0.3762  0.0000  0.0984  0.0026   0.8990
79  -1.9913 -0.5113  0.0000  0.0612  0.0182   0.9205
80  -2.6193 -0.4920  0.0000  0.0558  0.0257   0.9184
81  -3.4195 -0.6133  0.0000  0.0536  0.0177   0.9287
82   3.5971   3.4158  0.0000  1.0000  0.0000   0.0000
83   2.5976  2.5343  0.0000  1.0000  0.0000   0.0000
84   1.5995  1.5836  0.0000  1.0000  0.0000   0.0000
85   0.6033  0.6070  0.0000  1.0000  0.0000   0.0000
86  -0.2828 -0.1051  0.0000  0.4638  0.0007   0.5355
87  -1.0311  -0.3832  0.0000  0.0994  0.0000   0.9006
88  -2.3781  -0.8818  0.0000  0.0610  0.0001   0.9389
89  -3.3297 -1.2288  0.0000  0.0560  0.0001   0.9439
90  -4.2500 -1.5744  0.0000  0.0539  0.0000   0.9461
91   3.5487  3.3890  0.0000  1.0000  0.0000   0.0000
92   2.5490  2.5032  0.0000  1.0000  0.0000   0.0000
93   1.5500  1.5425  0.0000  1.0000  0.0000   0.0000
94   0.5541   0.5581  0.0000  1.0000  0.0000   0.0000
95  -0.2818 -0.1047  0.0000  0.4663  0.0003   0.5335
96  -1.0465 -0.3888  0.0000  0.0980  0.0000   0.9019
97  -2.4431  -0.9071  0.0000  0.0605  0.0000   0.9395
98  -3.4491  -1.2806  0.0000  0.0556  0.0000   0.9443
99  -4.4779 -1.6943  0.0000  0.0536  0.0000   0.9464
100  3.5000  3.3778  0.0000   1.0000   0.0000   0.0000
101   2.5000  2.5022  0.0000   1.0000   0.0000   0.0000
102  1.5000  1.5032  0.0000   1.0000   0.0000   0.0000
103  0.5000  0.5047  0.0000   0.9999   0.0001   0.0000
104  -0.2531  -0.0940  0.0000   0.5252   0.0009  0.4739
105  -1.0779  -0.4005  0.0000   0.0953   0.0003  0.9044
106  -2.4806  -0.9210  0.0000   0.0602   0.0000  0.9398
107  -3.4878  -1.2954  0.0000   0.0555   0.0000  0.9445
108  -4.5196  -1.7111  0.0000   0.0535   0.0000  0.9465

  FINAL VOLUME OF WATER =        17.0679
  FINAL VOLUME OF OIL  =        1.0120
                                                        93

-------
TOTAL ITERATIONS
TOTAL TIME STEPS
327
327
  ** RESULTS OF TRANSPORT *
  NODE-I-CONCENTRATIONS  IN PHASE-
     WATER  OIL  GAS    SOLID
SPECIES NUMBER = 1
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
0.4086E-03
0.1754E-02
0.2040E-01
0.1983E+00
0.1120E+01
0.1361 E+01
0.6695E+01
0.6556E+02
0.8976E+02
0.2163E-02
0.3079E-02
0.481 6E-01
0.5330E+00
0.451 8E+01
0.6243E+01
0.8732E+01
0.8939E+02
0.1135E+03
0.1739E-01
0.4171E-03
0.1478E+00
0.1675E+01
0.1647E+02
0.2443E+02
0.8065E+02
0.1654E+03
0.1632E+03
O.OOOOE+00
0.1397E+00
O.OOOOE+00
0.5898E+01
0.4600E+02
0.1347E+03
0.1781E+03
0.1778E+03
0.1789E+03
0.1762E+00
O.OOOOE+00
0.1514E+01
0.7630E+01
0.1478E+03
0.1812E+03
0.1774E+03
0.1780E+03
0.1778E+03
O.OOOOE+00
0.6265E+00
O.OOOOE+00
0.3002E+02
0.1797E+03
0.1768E+03
0.1781E+03
0.1778E+03
0.1779E+03
0.1306E+00
0.4056E+00
0.1074E+01
0.4291 E+02
0.1797E+03
0.1794E+03
0.1778E+03
0.1778E+03
0.1779E+03
0.2821 E+00
0.6066E-01
0.3535E+01
0.501 8E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1204E+05
0.3976E+05
0.8156E+05
0.8047E+05
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.2268E+05
0.6639E+05
0.8782E+05
0.8765E+05
0.8820E+05
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7287E+05
0.8932E+05
0.8744E+05
0.8776E+05
0.8767E+05
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.8859E+05
0.871 6E+05
0.8778E+05
0.8767E+05
0.8769E+05
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.2115E+05
0.8860E+05
0.8846E+05
0.8767E+05
0.8768E+05
0.8772E+05
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.3266E+00
0.1607E+01
0.1573E+02
0.2154E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1084E+01
0.1498E+01
0.2096E+01
0.2145E+02
0.2725E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.3954E+01
0.5863E+01
0.1936E+02
0.3970E+02
0.391 8E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1104E+02
0.3232E+02
0.4275E+02
0.4267E+02
0.4294E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.3547E+02
0.4348E+02
0.4257E+02
0.4272E+02
0.4268E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.431 3E+02
0.4243E+02
0.4273E+02
0.4268E+02
0.4269E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.431 3E+02
0.4307E+02
0.4268E+02
0.4268E+02
0.4271 E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
                                             94

-------
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
0.1657E+03
0.1777E+03
0.1768E+03
0.1778E+03
0.1779E+03
0.3233E+00
O.OOOOE+00
0.6116E+01
0.5632E+02
0.1356E+03
0.1653E+03
0.1745E+03
0.1763E+03
0.1767E+03
0.3442E+00
O.OOOOE+00
0.9120E+01
0.5930E+02
0.1187E+03
0.1454E+03
0.1639E+03
0.1710E+03
0.1705E+03
0.2995E+00
0.3389E-01
0.1154E+02
0.6081 E+02
0.1099E+03
0.1331E+03
0.1487E+03
0.1579E+03
0.1609E+03
0.2448E+00
0.1237E+00
0.1234E+02
0.61 21 E+02
0.1079E+03
0.1282E+03
0.1439E+03
0.1533E+03
0.1566E+03
0.8168E+05
0.8762E+05
0.871 8E+05
0.8767E+05
0.8769E+05
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.8147E+05
0.8601 E+05
0.8694E+05
0.8709E+05
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.3976E+02
0.4265E+02
0.4244E+02
0.4268E+02
0.4269E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.3255E+02
0.3966E+02
0.4187E+02
0.4232E+02
0.4240E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.2848E+02
0.3491 E+02
0.3934E+02
0.4105E+02
0.4092E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.2637E+02
0.3196E+02
0.3568E+02
0.3790E+02
0.3861 E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.2590E+02
0.3077E+02
0.3453E+02
0.3680E+02
0.3758E+02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
TOTAL MASS IN
  WATER    OIL     GAS     SOLID
0.5090E+03  0.8736E+05  0.4345E+03  O.OOOOE+00
SPECIES NUMBER = 2
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
0.2692E-03
0.1017E-02
0.5936E-03
0.9402E-02
0.3570E-02
0.2407E+00
0.4893E+00
0.5020E+00
0.4151E+00
0.3112E-02
0.2099E-03
0.2469E-01
O.OOOOE+00
0.2656E+00
O.OOOOE+00
0.1008E+00
0.6104E-01
0.1205E+00
O.OOOOE+00
0.2040E-01
O.OOOOE+00
0.391 6E+00
O.OOOOE+00
0.4633E+00
0.7895E+00
0.7501 E+00
0.751 1E+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.4633E+06
0.7895E+06
0.7501 E+06
0.751 1E+06
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.481 4E-03
0.9785E-03
0.1004E-02
0.8302E-03
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.531 2E-03
O.OOOOE+00
0.201 6E-03
0.1221E-03
0.241 OE-03
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.9266E-03
0.1579E-02
0.1500E-02
0.1502E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
                                                  95

-------
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
0.8927E-02
0.2759E-02
0.7571 E-01
O.OOOOE+00
0.5160E+00
0.7644E+00
0.7442E+00
0.7443E+00
0.7438E+00
0.5189E-02
0.1095E-01
0.4728E-01
0.3296E-01
0.7583E+00
0.7428E+00
0.7445E+00
0.7442E+00
0.7443E+00
0.2833E-02
0.1747E-01
0.1 151 E-01
0.1522E+00
0.7435E+00
0.7448E+00
0.7442E+00
0.7443E+00
0.7443E+00
0.5655E-02
0.1 251 E-01
0.2329E-01
0.7950E-01
0.7435E+00
0.7436E+00
0.7443E+00
0.7443E+00
0.7443E+00
0.6608E-02
0.1068E-01
0.2964E-01
0.3954E-01
0.6200E+00
0.7444E+00
0.7448E+00
0.7443E+00
0.7443E+00
0.7669E-02
0.8630E-02
0.3527E-01
0.8822E-02
0.1254E+00
0.7548E+00
0.7459E+00
0.7450E+00
0.7449E+00
0.4267E-02
0.1618E-01
O.OOOOE+00
0.1639E+00
0.1307E+00
0.6823E+00
0.7535E+00
0.7449E+00
0.7470E+00
0.8855E-02
0.7154E-02
0.3209E-01
0.1682E-01
0.8270E+00
0.7555E+00
0.7452E+00
0.7471 E+00
0.7459E+00
0.7023E-02
0.1 071 E-01
0.2120E-01
0.3906E-01
0.7671 E+00
0.7267E+00
0.7352E+00
0.741 9E+00
0.7479E+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.5160E+06
0.7644E+06
0.7442E+06
0.7443E+06
0.7438E+06
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7583E+06
0.7428E+06
0.7445E+06
0.7442E+06
0.7443E+06
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7435E+06
0.7448E+06
0.7442E+06
0.7443E+06
0.7443E+06
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7950E+05
0.7435E+06
0.7436E+06
0.7443E+06
0.7443E+06
0.7443E+06
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.6200E+06
0.7444E+06
0.7448E+06
0.7443E+06
0.7443E+06
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.7548E+06
0.7459E+06
0.7450E+06
0.7449E+06
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1032E-02
0.1529E-02
0.1488E-02
0.1489E-02
0.1488E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1517E-02
0.1486E-02
0.1489E-02
0.1488E-02
0.1489E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1487E-02
0.1490E-02
0.1488E-02
0.1489E-02
0.1489E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1487E-02
0.1487E-02
0.1489E-02
0.1489E-02
0.1489E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1240E-02
0.1489E-02
0.1490E-02
0.1489E-02
0.1489E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.2509E-03
0.1510E-02
0.1492E-02
0.1490E-02
0.1490E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.261 4E-03
0.1365E-02
0.1507E-02
0.1490E-02
0.1494E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1654E-02
0.1511E-02
0.1490E-02
0.1494E-02
0.1492E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
0.1534E-02
0.1453E-02
0.1470E-02
0.1484E-02
0.1496E-02
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
O.OOOOE+00
96

-------
  TOTAL MASS IN

  WATER   OIL     GAS     SOLID
0.2229E+01 0.7511E+06  0.1637E-01  O.OOOOE+00

END OF SIMULATION RUN
                                                97

-------
APPENDIX D. FLOW CHART FOR PROGRAM
                98

-------

-------
                    APPENDIX E: DESCRIPTION OF SUBROUTINES

ARNG       Arranges active nodes for ASD options in ascending order

BMASS      Computes volume of fluids in various phases from saturations for Cartesian domain

BONEL      Sorts out elements on the boundary of the domain

CALS        Returns the local and the global node number of nodes given the element number
             and the side for assigning type-3 transport boundary conditions

CONN       Gives numbers of connected nodes for any node in the domain

DENU       Updates  fluid densities from fluid composition

DLST2       Computes element stiffness matrix for transport problem

DMASS      Computes the mass for each species in different phases

DOB V       Computes type-3 boundary conditions for transport problem for an element side
             from the assigned schedule

DOBW       Computes fluid fluxes for type-2 flow boundary conditions for an element side
             from the assigned schedule

DODIF       Imposes type-1 boundary conditions for transport problem

DVELO      Computes fluid velocities for water, oil and gas phases at nodes in the domain

EJACK       Evaluates element stiffness matrix for flow problem

ELST4       Assigns area for Cartesian problem or volume for radial problem to each node in
             the domain

FLUX        Calls DOBW for an element with type-2 flow boundary conditions and modifies the
             element stiffness matrix accordingly with the information returned from DOBW

GASX       Assigns the prescribed minimum gas pressure to nodes with total liquid
             saturation =  1

GATE1       Evaluates element properties required in DLST2

JMATER     Computes nodal properties required in the element stiffness matrix for flow
problem

                                         100

-------
MASSCH    Computes mass injection rates for components in the transport module

MJACK      Assigns nodal properties computed in JMATER for each node of the element to the
             relevant arrays and sends to EJACK

MTIPH      Computes interphase mass transfer terms

RDGAS      Reads input data for flow module

RDTR       Reads input data for transport module

REDBC      Imposes boundary conditions for the Newton-Raphson scheme

REDNSY    Solves the global assembled equations

RMASS      Computes the volume of fluids in various phases in a radial domain

RWRO       Modifies load vector in the fluid flow module to account for interphase transfer

SODNW     Assigns values of type-1 flow boundary conditions to nodes from the corresponding
             schedules

TRANSP     Main control subroutine for the transport module

UPFAM      Computes nonequilibrium interphase partition coefficients

UPWIND2    Computes the upstream weighting for the flow module

WRTR       Prints output for the transport module

XYZ         Computes nodal coordinates and node numbers for each element in the domain
                                         101

-------
               APPENDEX F DEFINITIONS OF IMPORTANT VARIABLES
ABSA       Absolute convergence limit for all liquid heads

ABSC       Absolute convergence limit for concentration of species in water phase under
             nonequilibrium mass transfer conditions

AKWX      Conductivity of water in x-direction

AKWY      Conductivity of water in y-direction

AKOX       Conductivity of oil in x-direction

AKOY       Conductivity of oil in y-direction

AKAX       Conductivity of gas in x-direction

AKAY       Conductivity of gas in y-direction

ALM        Gas compressibility

A           Global stiffness matrix

B           Global load vector in the flow module

AMASS      Accumulated mass for species in various phases

ACN        Boundary condition schedules in the transport module

ACDF       Node numbers connected to nodes in the domain (connectivity matrix)

BRN        Temporary storage for material properties in the flow module

CFn         Temporary storage for element material properties required in the evaluation of the
             element stiffness matrix in subroutine EJACK (n=l, 2, 3, 4, 41, 42, 44, 45, N41,
             N42, N44, N45)

CHU        Array to store fluid pressure heads at which derivatives in the Newton-Raphson
             scheme are evaluated

CHA        Array to store fluid pressure heads at start of the time step

CHAF       Array with current updated fluid pressure heads

                                          102

-------
CHAI        Previous iteration fluid pressure heads




CON        Species concentrations at the start of the time step




CONF       Species concentrations at the end of the time step




DELT        Time increment




DTMX       Maximum allowable time increment




DTMIN      Minimum allowable time increment




DCAY       First order decay coefficient for species in various phases




DENI        Density of the inert organic component




DIFW        Diffusion coefficient of species in water




DIFO        Diffusion coefficient of species in oil




DIFA        Diffusion coefficient of species in air




DLIQ        Liquid density of pure species




UDENW     Updated density of water phase




DENO       Updated density of oil phase




DENA       Updated density of gas phase




DXX        Dispersion coefficients in the x-direction




DYY        Dispersion coefficients in the y-direction




ECORD(i,l)  x-coordinates of nodes i in the domain




ECORD(i,2)  y-coordinates of nodes i in the domain




EPART      Nonequilibrium partition rate coefficients




ESM        Temporary storage during element stiffness matrix computations




EF           Temporary storage during element stiffness matrix computations






                                          103

-------
ESN         Temporary storage during element stiffness matrix computations

ESN1        Temporary storage during element stiffness matrix computations

ESMN       Temporary storage during element stiffness matrix computations

ESMN1      Temporary storage during element stiffness matrix computations

ESMNX     Temporary storage during element stiffness matrix computations

GAMA      Equilibrium partition coefficients between phases (see App. A)

HA          Nodal gas pressure at the start of the time step

HAF         Updated gas pressure at nodes

IRUN        Index for execution control (see App. A)

IRES        Restart option index (see App. A)

ICOM        Index describing gas compressibility (see App. A)

ICON        Index for convergence 1  = no convergence and the iteration is repeated;0=solution
             has converged

ICN         Number of flow BC schedules

IDIM        Index to specify the units of linear dimension (see App. A)

IPROP       Material property to be assigned to the node

IRAD        Index to describe 2-D flow (see App. A)

ITRN        Index to include transport module (see App. A)

IDBC        Array to store number of element and sides subjected to type-3 BC transport
             module (see App. A)

IFEMC      Index to activate ASD (see App. A)

IMCOR      Index to activate mass injection scheme (see App. A)

IMPE        Serial number of an active element in ASD formulation
                                          104

-------
JINF
Index for terminating the current run based on cumulative fluid intake
KEQUL      Index to activate kinematic mass transport (see App. A)
KUNIT      Index to specify mass units (see App. A)

KTUDAT    Index to specify updating of mass transfer rates (see App. A)

KCN         Total number of schedules describing transport boundary conditions

KIC         Option for initial condition of transport (see App. A)

KSCH       Schedule for transport boundary conditions

MCNT       Maximum number of iterations allowed

MXTST      Maximum number of time steps allowed

NCM        Node numbers in an element

NPH         Number of phases (see App.  A)

NNBi        Nodes and schedules for water (i=l), oil (i=2) and gas (i=3) phases for type-1 flow
             boundary  conditions (see App. A)

NNSp        Elements  and side numbers for water (p=W), oil (p=O) and gas (p=A) phases
             subject to type-2 flow boundary conditions (see App.  A)

NTSP        Total number of noninert species

NCON       Total number ot type-1 nodes for transport

NEBC       Total number of type-3 elements for transport

NMAT       Total number of material types

NSPB        Number of node with type-1 transport boundary condition

PROP        Material properties (see App. A)

PRON       Temporary storage of element material  properties in flow module

PR          Temporary storage of nodal properties in flow module

PHI         Porosity

                                         105

-------
vxx
x-direction flow velocities at nodes
VYY        y-direction flow velocities at nodes

QFX        Subset of x-direction flow velocities at nodes in an element

QFY        Subset of y-direction flow velocities at nodes in an element

RHOW      Specific gravity of oil phase

RKNA       Relative conductivity of phase at nodes

RHP( 1,1)    Density of pure water

RHP(1,2)    Density of inert oil component

RHP(1,3)    Density of uncontaminated air

RW(i,l)      Interphase transfer rate to water phase at node i

RW(i,2)      Interphase transfer rate to oil phase at node i

RW(i,3)      Interphase transfer rate to gas phase at node i

SAW        Updated saturation of water phase

SAO        Updated saturation of oil phase

SAP         Subset of liquid phase saturations

S AWO       Water phase saturation at the start of current time step

S AOO       Oil phase saturation at the start of current time step

SWMN1     Historical minimum saturation at a given location since changing from a two phase
             air-water to a three phase air-oil-water system

SOT         Saturation of trapped oil

TINF        Cumulative total intake (see App.  A)

TIME        Simulation time

TMAX       Maximum simulation time

                                           106

-------
TMXG       Tolerance for deciding if gas phase is at steady state (see App. A)




TH          Time weighting factor (see App. A)




TFACT      Increment factor for time-step (see App. A)




TOLL        ASD tolerance for liquid pressures (see App. A)




TOL2        ASD tolerance for liquid saturations




TPNT        Time interval for printout of results (see App. A)




TRANS(l)   Longitudinal dispersivity




TRANS(2)   Transverse dispersivity




TTRN        Time at which interphase mass transfer updating starts




URW        Interphase mass transfer rate between oil-water phases




VRA        Interphase mass transfer rate between air-oil/water phases




VRS         Interphase mass transfer rate between water-solid phases




VISAW      Ratio of air to water viscosity




VOL        Area (Cartesian) or volume (radial) of the domain assigned to each node




WFAC       Mass injection rates at nodes




WFNn       Temporary storage for element properties in flow module (n=41, 42)




XCO        x-coordinates of nodes on the x-axis starting from the lower left corner




ZCO        y-coordinates of nodes on the y-axis starting from the lower left corner
                                          107

-------
                               O)
                               CD
Figure 4,L Local side mmibering of quadrilateral
                                               6ft    72
                      (93
Figure 4,2, Finite e-lcaiaat mesh with numbers, m parcntheses indicating
          global element numbefe and dChecs referda^ to global
          tiodi!

-------
                           I
                             tHM>-
                            u
                               C_Q
                                                 EthinUwi
                              JCtti.fi
                              «0-.D
                                  •  • ti PI i ii ij 111 ii im LI 11 L ii  111 M i • ii i ii i i
                                 tutfir    dJi>    HUM    (Ln    flJb    1J»
Figure 5.1  Witter and total liqwid saturadofl v» depth At end of Stages L, 2 md 3 fa Eseample 1-

-------
                       Time (DA?*)
5-2
*
   ftf-
     M
     ft 00
   for
icyl^e Bid tol««

     Uld
                                      t outflow lx«ad»ry fa Eumpk 1

                             analyses. Rate OOcfficknts in d'1.

-------
                                                          11-fl
                              K(m)
Figure 5-3  Qit Mturs-tio-n conioura it the end of Stife 2 far Example 2,

-------
 M
                                                            JLfr
  fed    1.D
     5.4  Water phise
                                 [j m'1) at tlie end of Stage 2 for Example 2-
  M   «   U   4.0   -M!    S.O   «,»   3J   M   M   Ifeff  TtJt_[
                  uimii	mil	•MMMMMMMirM"-—	I  ' • ^^^^wl    t ji  MBBBBBBBBPFJ--^ M    1 ^,|
B..II
TJ>  -
  6d    1.0   M   SA   **
Fi|ure 5.5 Ga* phi*e «m«ntrttioii contoun (j raf*) at the end «f Sla^e 2 tor
                                                                         2.

-------
                                           15J   **    rtV
   SCO
      «     1W
                               R{cm)
       5.6  Oil saturation ooatouts at tie cad of St»gE 2 for Exunplfi 3-
   3KU3
      1ftfl    IWfl   1HW
 £
 u
     DJ
                   llftfl
                                 R(crn)
                                 1 - " - 1 - 1
                                                                WO.D-
Fi§ur« 5.1  Water BaluratiMi contour* At the end of Stage 3 for Exainpl* I.

-------
' ' i ' fc "• j '  * '  J i*  J '  '
4 SOE*006 *^"t i i
         0,00
    5.&  PER. mass in system v§ time during Stage 3 of Esaniple 3-

-------