Immisicble two-phase flow drainage in a subsurface reservoir - example 1

Load and plot the petrophysical (permeability and porosity) data

We considere a simple one layer rectangular domain. In this case we explicitly set the domain lengths along each space dimension as follows

Next, we set the number of cells along each space direction and the total number of grid cells. Next, we add those variables into the Grid structure

To fullfill the need of internal compuattional tasks we add the uniform grid spacings, hx, hy, and hz, along the respective directions to the Grid structure

Additionally, we need to setup other few variables for the Grid structure such as, V (the uniform cells volume), and compr (the total compressibility). Notice that the compressibility data is a 3D shape array whose size conforms to the grid dimensions. The MATLAB command ones(Nx,Ny,Nz) constructs such an array with a unit values for all cells, and thus by multiplying it with a constant values we obtain a uniform distribution of the compressibility data over the computational domain.

Now, we load porosity data from a previously saved mat file 'spe_layer1_phi.mat' which is included in the data subfolder of the included files. porosity data is stored in phi_layer variable (because we can store more than one array in a mat file).

Now, we will check if the loaded data is within acceptable bounds. A situation which occur many often in practice. First, look if there exists porosity values below or close to 0+. This can be simply done with the following command:

indicating that the porosity is below the smallest floating point precision number eps (effective zero on the target computer) at 68 cells. To correct for this misleading data we set these values to a prescribed minimum porosity threshold (1e-3 herein). Otherwise, the saturation solver will fail since dividing by zero porosity will result in an overflow.

The Grid structure is promoted with a new variable por conforming to the Grid shape.

Finally, permeability data is loaded from the mat file 'spe_layer1_perm.mat' and variable perm_layer.

Tha data units is converted into SI unit of the permeability (m^2) because it is given in Dary units in the SPE10 dataset. Again, the Grid structure is promoted with a new variable K conforming to the Grid shape but with 3 components since the permeability data assumes an anisotropic diagonal tensor.

Now, we're ready the plot the porosity and permeability data to see how they look like. First, we will plot the porosity data in linear scale. Next, we plot permeability along X direction in log10 logarithmic scale to depict its extent over several orders of magnitude.

To do this, a cartesian grid G data structure is setup from the domain extent ans spacings. This will help to do the plotting commands more efficiently.

Start a new MATLAB figure and specify the colormap corrsponding the color palette to use and the number of colors in this palette.

Plot the porosity data into a subplot in the same window but at le left hand side. The PlotCellData function plots a map of the variable Grid.por using information in the structure variable G. The 'EdgeColor' option is set to 'none' to hide the display o the gridlines of the underlying grid. Since the grid is so dense in our modl, letting this option active the resulting plot will be obsured.

next, this plot x and y axis are set as dependent and we add a color bar and a title.

The same sequence of commands are repeated for the right plot to obtain a map of Kx in log10 scale.

subplot(1,2,2);

PlotCellData(G, reshape(log10(Grid.K(1,:,:,:)),N,1), 'EdgeColor','none');

axis tight equal; colorbar, title('log_{10}(K_x) permeability map [m^2]');

Clear phi_layer variable from memory space since it's not needed further. This is more relevant for large scale applications or on computers with very modest ressources.

Now, we check if the log10-permeability and porosity fields ara correlated somehow. First, let's start with the first component of the diagonal permability tensor, Kx:

Off-diagonals of the correlation matrix are high positive values indicating a good positive correlation between the two variables. Let's now do the same thing for the y and z components of the permeability to conclude on their correlation with the porosity.

Setup auxilliary simulation data

We further add a five-spot well pattern composed from 1 injector and 4 producers. Production wells are placed at the domain corners while the injection well is located at the center. A quarter of the volumetric injection flow rate, Qi, is produced at each of the four production wells.

Notice that (by convention), flow rates have to be positive for injection wells and negative for production wells. Qw is a distributed 1xN array which holds distributed flow rates. It has zero values everywhere except at injection/production cells.

Next, we initialize two fluid objects oil and brine in which fluid properties data corresponding to a more dense and viscous Oil phase and less viscous and dense Saltwater phase, respectively.

Initialize fluid global pressure in the computational domain to a constant value of 100 bars and the initial brine saturation to its residual value (i.e. brine.Sr)

Define time stepping parameters including total simulation time, T, the number of time steps, nt, and the global time step, dt.

For the nonlinear solvers we need to pass an optional structure with options holding the absolute tolerance, tol, and the maximum number of allowed Newton-Raphson iterations, maxiter. The latter corresponds to the inner iteration count since the most outer loop corresponds to time stepping. The last option, min_dt, is the minimum allowed time step size for the saturation solver. This is needed because, as explained later in more detail, we use an adaptive time stepping algorithm which cancels out if the time step size becomes smaller than this prescribed threshold.

Run the time simulation loop - Sequential implicit solution

Now, we will build the time loop involving a sequential solution of the governing global pressure equation and the nonlinear non-wetting phase saturation equation formulated in fractional form. The pressure and saturation solvers are both implicit implying an unconditional stability regardlless of the magnitude of the selected time step.

The pressure equation, neglecting the gravitational and capillary pressure terms is given as:

where is the total compressibility, is the rock porosity, is the fluid pressure, is the total mobility given as a sum of oil and water mobilities, respectively. A fluid phase mobility is defined as the ratio of its relative permeability to its viscosity . is the water injection/withdawal mass flow rate [M/T] and is the density of the water phase [M/L^3]. Threfore is the volumetric flow rate.

The saturation equation is a nonlinear diffusion-convection like equation for the water phase given as:

where is the water phase saturation, is the saturation-dependent fractional flow function for the water phase, is the total velocity (i.e. the sum of fluid velocities of each phase).

Herein, the pressure solver is linear but requires successive updating of the total mobility of both fluid phases which is the task of TwoPhasePressure function. The latter wraps the Two-Point Flux Approximation solver (TPFA function) such that is iteratively updated from spatial distribution of water saturation calculated in the previous time step.

The saturation solver uses an heuristic adaptive time steeping procedure combined with an inner Newton-Raphson (NR) iteration loop. During the inner NR iterations, the sparse system Jacobian matrix is evaluated analytically to accelerate the algorithm convergence. If the NR iteration fails the current time step is subdivided (i.e. the number of local time steps is multiplied by 2) and the whole time step is repeated again from the beginning. The latter mechanism is called back-stepping. The local NR failure is due to either (i) an insuffisant or small number of allowed maximal iterations to achieve the required accuracy. Hence, when setting the NR options the user should keep a good balance between the maxiter and tol options; (ii) second, to a high or very difficult level of prescribed accuracy which cannot be achived for badly preconditioned problems. So, the absolute tolerance should be kept to reasonable values (in practice values between 1e-3 and 1e-4 are suggested); (iii) and finally, when the selected global time step results in large local pressure jumps, the same time step is generally inappropriate for the saturation solver implying a smaller time step advancement of the saturation solver to absorp the sudden numerical shock caused by the large pressure jump.

In general, this adaptive time stepping strategy is effective in practice. Very often the algorithm starts with smaller saturation time steps at the begining of the simulation when initial conditions are far from the expected solution. This sub-time step will gradually increase in the next time levels if there is no abrupt changes in the boundary conditions (i.e. wells pumping or production schedules) or materials and fluids properties.

Let's run this two-phase flow oil drainage simulation

From the convergence output we notice that this sequential solution algorithm went very smoothly. There was a difficulty to reach convergence in the first time step because it seems that 2 days time step was very large for the saturation solver. Hence, the initial time step was automatically subdivided 11 times resulting to an effective saturation time step close to 84 seconds as shown below from the fields of the returned output structure out for the first time period. The dt field returns the selected time step by the saturation solver during this period.

Now, we plot a summary of the saturation solver performance during the simulation. In the left side of the figure we plot the evolution of the internally selected time step versus time, while in the right figure the evolution of the total number of inner Newton-Raphson iterations versus time is plotted. Both figure are best depicted using a semilogy axis.

figure;

subplot(1,2,1);

for t=1:nt, dt_history(t) = out(t).dt; end

semilogy((1:nt)*dt/day, dt_history/day, '-ro');

title('Evolution of the saturation time-step'); xlabel('Time [days]'); ylabel('Time step [days]'); grid on;

subplot(1,2,2);

for t=1:nt, nr_history(t) = out(t).nr_iterations; end

semilogy((1:nt)*dt/day, nr_history, '-b*');

title('Newton-Raphson iterations'); xlabel('Time [days]'); ylabel('Total NR iterations'); grid on;

We observe thet the internal saturation time step gradually increases to 1.5 hour, 6 hours, 12 hours and finally to 2 days starting from 22 days after the simulation start. So, afterwards the pressure and saturation time steps becomes equal meaning that the saturation extent of the plume becomes significant enough to amortize the pressure jumps immediately.

Post-processing of the numerical solution - mapping mode

Since the main simulated variables are pressure and water saturation we will plot them directly using the supplied functions which make use of MATLAB's excellent built-in graphics capabilities without writing the results into intermediate text or binary files in a required format and loading them back from an external third-party post-processing software. So, with very few commands we get directly the results ready for quality publication freeing the modeller from all intermediate tasks and hence speeding the data-processing-simulation-post-processing workflow significantly.

We will use here the routine PlotCellData to produce the figures in mapping mode where each cell in the grid is plotted with a constant color corresponding to a linearly interpolated value between the upper and lower colors in the used color palette. This palette fits the range between the maximal and minimal values detected in the data array.

The figure command without any argument opens a new MATLAB figure (or axis in MATLAB's terminology) so that the previously opened figures are preserved and all following drawing commands will be placed in this new figure. The colormap command is optionally called before the drawing commands to choose among the available color palettes, here we choose the jet palette with 256 colors distributed between blue and red.

Next, we plot the pressure after conversion from Pascal to bar units. We call the subplot command to specify that we want to make a plot matrix with 1 row and 2 columns in the target figure. The last argument in the subplot command is the active subplot to which the next plotting commands are directed. The PlotCellData function takes additional property-value pairs (à la MATLAB). These options are exactly those supported by the low level patch command to control the internal details of the plot. To study all these options type help patch command or go to patch window in the online help system. Here the 'EdgeColor' property specifies the color of the underlying grid edges, by setting this property to 'none' they will be hidden in the current plot. In this example the grid is very dense. So, drawing all grid edges will obscure the reulting plots but in other applications where the grid is coarse this might not be necessary.

By default MATLAB plots fit to the rectangular window inside which they're displayed. A good idea here is to call the axis command to make x and y axis scales dependent using equal option. Likewise, tight wil readjust the plot to fit at best to actual window after axis scaling. Finally, we decorate the plot with a colorbar and a custom title.

We follow the same previous steps to plot the water saturation in our computational domain.

% plot brine saturation on the right side of current figure

subplot(1,2,2);

PlotCellData(G, S(:), 'EdgeColor','none');

axis tight equal; colorbar, title('Brine saturation [-]');

Post-processing of the numerical solution - contouring mode

Presenting the results in mapping mode is most natural when working with cell-centered discretization methods. It has also the advantage to appreciate the dispersion of the saturation plume as seen in the right plot of the above figure. But, it has other limitations. For instance, the gradient of the plotted variable cannot be easily appreciated in some areas, the presentations graphics quality may deteriorate when saving the figure into raster file formats, and this presentation mode is not suitable for explanation to non-specialists such as managers and decision makers.

Another presentation mode which is quite popular in all fields of geosciences in contouring. The target variable iso-values are traced by a general purpose algorithm to contour the data. Such presentation mode is mainly supported by basic MATLAB routines such as contour, contourc, contourf and contour3. We will use PlotContourData function which is a simple wrapper to contourf or contour to plot filled or simple contour plots, respectively. This function arguments are similar to those of PlotCellData function except for the optional arguments are property-value pairs which are compatible with those of the contourf or contour commands.

So, make a figure made from two contour plots of water pressure and saturation. We plot pressure contour lines between 20 and 140 bars with an interval of 20 bars. We select saturation contour lines between 0 and 0.7 with 5% interval.

figure;

opts.fill = 'off'; % change this option to 'on' for contour filled plots

% plot water pressure on left side

subplot(1,2,1);

PlotContourData(G, P(:)/bar,opts,'LevelList',20:20:140);

colormap(jet(16));

axis tight equal; colorbar, title('Pressure [bar)]');

% plot water saturation on right side

subplot(1,2,2);

PlotContourData(G, S(:),opts,'LevelList',0:0.05:0.7);

axis tight equal; colorbar, title('Brine saturation [-]');

To better understand the spatial distribution of brine saturation we also plot its probability distribution reflecting the saturation volumetrics inside the reservoir. To do that we call the histogram command for variable with 7 bins as a second argument and 'probability' value for the 'Normalization' option. This means that the height of each bar is the relative number of cell-centered saturation items (number of saturation values in bin / total number of grid cells). The sum of the bar heights is obviously equal to 1.

figure;

histogram(S,7,'Normalization','probability');

title('Saturation probability distribution');

Around 35% of the reservoir volume is filled with water saturations between 0.4 and 0.5 while around 75% of the reservoir volume is filled with water saturations between 0.3 and 0.6 which are clearly the dominant modes of saturations after 50 days of oil drainage by water injection.

Conclusions

In merely a little bit more than 100 lines of MATLAB commands (including blank and comment lines) we were able in this tutorial to build a script performing the following tasks:

• Read petrophysical data from disk files, plot, and analyse them.
• Setup auxilliary simulation data.
• Run an immiscible two-phase flow drainage simultion involving an implicit sequential simulation procedure, analyse and plot the convergence behavior of the adaptive time stepping algorithm.
• Post-process the numerical solution using mapping and contouring modes, then analyse the spatial distribution of water saturation after 50 days of oil drainage.

All these steps take less than 1 minute CPU time in a commodity PC. With a scripting technology the time traditionally taken in file and data transfer between pre-processing, simulation, and post-processing tools as in existing workflows is cutted to its minimum. Hence, the modeller can spend more time to gain in problem understanding rather than loosing a significant project timeframe in useless data manipulation tasks.