The actual program has been coded for a two-dimensional system, but can easily
be extended to three dimensions, as the equations have shown.
All equations are solved implicitly for the pore pressure *P*, the gas
saturation *S*_{g} and the hydrate saturation *S*_{h} at each grid point with the
Newton-Raphson method using a Cartesian grid system.
A flow chart of the code can be seen in Figure 4.

The time is denoted by *n*, the number of Newton iterations by *it*, and
the time steps by .
Each new is calculated
automatically based on the
following relationship Aziz (1995):

(30) |

The desired change of either pore pressure *P* or saturations *S*_{g} and *S*_{h}
over is given by , the actual change by and
is a tuning factor with values between zero and unity. Usually,
has the same value for all grid blocks, but different values for the
different variables (in this case, pressure *P*, gas and hydrate saturation
*S*_{g} and *S*_{h}). So a different is calculated for each
variable and the minimum is used for the subsequent time step.

Figure 4

Initially, the system is characterized by a given
pore pressure, temperature, gas and hydrate saturation. Perturbation of the
pressure at time *t*=0 or an initial system, which is not in equilibrium,
will initialize fluid-flow through the system. The boundary conditions
depend on the initial system, and assume either constant pressure
boundary at *z*=0 and no-flow boundaries everywhere else,
or no-flow boundaries everywhere.

The code is still in a test phase. As of now, it is working for only two-phase flow and such small hydrate saturation that neither hydrate growth nor dissociation can occur. It is still unstable and non-converging for hydrate saturations bigger than 0.01. Additional work is required to correct this problem.

I conducted a simple test on the functionality of the code for at least two-phase fluid flow with the third, solid hydrate phase being approximately zero. The initial system consisted of a no-flow boundary setup of water over gas with a negligible small hydrate saturation everywhere. The initial pore pressure was comparable to hydrostatic pressure, and the temperature was chosen so that neither hydrate growth nor decomposition could be expected for the pore pressure range. The density difference of gas and water, the pressure gradient and the effect of gravity result in a fluid flow stimulation in which the gas migrates upward and the water travels down.

The grid system is a 10x10 system, each grid block having dimensions of
10 m x 10 m. The top of the grid system is at about 1000 m depth. The grid
blocks within the first 40 m are filled with 100% water, the last 60 m
are filled with 100% gas. To avoid any instabilities resulting from
a vanishing third
flow equation (hydrate equation), a constant hydrate
saturation of 0.001 was assumed throughout the system. The temperature
was chosen to be 10 degrees, and the initial pore pressure followed the
relation *P*_{ij}= 1000 kg/cm * 9.81 m/*s ^{2}* *

Figure 5 shows the gas saturation results of the fluid-flow initialized by the instable gas-water system. 100 % gas saturation is displayed as white, 0 % as black. At time 0, the water is above the gas and fluid flow is stimulated. After 1 day, gas and water have started to migrate depending on low and high permeabilities in the different grid blocks. Following the permeability pattern, the two fluids begin to mix. After 155 days, the two phases have nearly equilibrized and finally they reach their equilibrium state with gas above water. As soon as this equilibrium state is reached, there is no more flow through the system.

Figure 5

The test shows that the code appears to work correctly for at least a two-phase fluid flow. More work is still required to improve upon the incorporation of a third-phase with more than 0.001 saturation everywhere.

11/12/1997