Eclipse Run Problems And Errors

Here is discussed ECLIPSE run problems and errors; what the reasons and possible solutions. Most of the advice and suggestions apply equally to all the simulators.
Note: To know about Eclipse Severity Levels read the article

Convergence in Eclipse

Eclipse starts to solve the equation(Solve = find the Pressure(P) and Saturation(S)) base on P and S value in the previous time step. Eclipse determines the length of time step by an internal algorithm based on how difficult the equation is solved in previous time steps.
The solution mechanism in eclipse is "Try and Error", the first guess is the values of P and S in previous time step. Therefore, the smaller the time step, the closer the initial guess to final value (smaller time step= smaller change in P and S)
If Eclipse fail to solve the equation for the chosen time step, it will give you a problem and chops the time step. It works sometime, but not always.
If time step chopping doesn't work, Eclipse gives you the error "Time step is too small to chop".

See Background to timesteps and iterations

Non-Linear Convergence Criteria

Non-linear iteration

The equations that the simulators are trying to solve are non-linear. By non-linear we mean that for instance doubling the tubing-head pressure of a water injector will not usually double the amount of water injected, and doubling the oil saturation in a grid block will not usually double the oil mobility in that grid block.

The simulators use an iterative process based on Newton's method to solve these non-linear equations:

We linearise the equations
We solve the linear equations
We check if this linear solution gives us a good enough non-linear solution.
If it does then we move to the next timestep. If not we calculate the change needed to improve the solution, then go back to step 1.

There is a limit to the number of non-linear iterations that the simulators will try before giving up and trying with a smaller timestep. This limit depends on the simulator and on the solution method, and is set using either the TUNING or the CVCRIT keyword.
You can check on how well the model is converging to a solution by looking at the number of non-linear iterations for each timestep, as described here:

1 non-linear per timestep means the step was very easy to converge
2 to 3 non-linears per timestep means the step was easy to converge
4 to 9 non-linears per timestep shows an increasingly difficult problem
> 10 non-linears per timestep can mean a problem with the model

Definition of convergence

Formally, we need a test to decide when we can stop the iterations and carry on to the next timestep, i.e. to decide when the solution that we have is 'good enough'. A test of convergence can be based on either a small residual or a small change in solution (or both).
The convergence criteria are different in ECLIPSE 100 and ECLIPSE 300.
The residual can be thought of as a measure of how close we are to solving the non-linear problem and is used in ECLIPSE 100 for the primary test. ECLIPSE 300 uses the change of solution as its basic convergence test. In the case of Dual Porosity ECLIPSE 100 runs, the TUNINGDP keyword will accept iterations as converged if either the residual or the solution change is small enough. The TUNINGDP keyword is sometimes helpful for high throughput cases.

Solution Change Check

The variables used by ECLIPSE 300 are pressure and molar densities. ECLIPSE 300 calculates the maximum change over all cells for these variables and compares these maximum changes to two limits, one for pressure (SCONVP) and one for saturation (SCONVS).

SCONVP Set by 1st data item in CVCRIT default 0.1 atm
SCONVS Set by 7th data item in CVCRIT default 0.01

The pressure change is compared directly with SCONVP, while an approximation using field-wide properties converts molar density changes to effective saturation changes. These effective saturation changes are compared with SCONVS. The maximum solutions change over all cells and all solution variables is converted to an "actual-over-target" (aot) which is a ratio of the current maximum divided by the appropriate limit. When the value of aot falls below 1, then the timestep is converged.
For example, suppose that the largest absolute pressure change in any grid block during the last non-linear iteration was 0.25 atm, with a limit on pressure change SCONVP = 0.1 atm. The ratio of the two number is 2.5, meaning that the pressure change was 2.5 time bigger than we would accept. Suppose the limit on effective saturation change was 0.01 and the actual largest absolute value in any cell was 0.02, which means the actual change was twice what we would accept. The value aot is the maximum of 2.5 and 2, so aot=2.5.
If the 8th data item in the DEBUG3 keyword is set greater than 0, then ECLIPSE 300 will output a line of debug that shows the convergence behaviour something like:

NLStep= 0 lin= 23 aot= 97.21 Rmax= .7162E-01 Rsum= .9919E-05 egain=-.1000E+01
NLStep= 1 lin= 19 aot= 17.55 Rmax= .1762E+00 Rsum= .1329E-06 egain=-.1000E+01
NLStep= 2 lin= 21 aot= 2.94 Rmax= .1421E-01 Rsum= .6285E-06 egain=-.1000E+01
NLStep= 3 lin= 12 aot= .77 Rmax= .7252E-02 Rsum= .7476E-08 egain=-.1000E+01
Rep ; 8901.0 1.00 8.7838 .19498 1.4E05 32884. 1.2E06 4843.6 .00000 1.3E06 4 2%
NLStep= 0 lin= 27 aot= 137.50 Rmax= .7587E-01 Rsum= .2011E-04 egain= .1805E+00
NLStep= 1 lin= 26 aot= 79.23 Rmax= .7589E-01 Rsum= .1743E-04 egain= .1676E+00
NLStep= 2 lin= 26 aot= 76.18 Rmax= .7279E-01 Rsum= .1676E-04 egain= .2621E+00
NLStep= 3 lin= 24 aot= 9.30 Rmax= .9062E-01 Rsum= .8301E-06 egain=-.1000E+01
NLStep= 4 lin= 24 aot= 9.00 Rmax= .8764E-01 Rsum= .8028E-06 egain=-.1000E+01
NLStep= 5 lin= 13 aot= .08 Rmax= .9509E-01 Rsum= .9853E-09 egain=-.1000E+01
MIF ; 8903.0 2.00 8.7860 .19501 1.4E05 32880. 1.2E06 4844.2 .00000 1.3E06 6 2%

This shows two time steps:
The first 4 lines above show iterations 0 to 3. Non-linear iterations start with iteration 0, which is a first guess at the new solution. The aot for iteration 0 is 97.21 so the non-linears are not converged.
Iteration 3 has an aot of 0.77, which is less than 1. The largest absolute change is now less than the target, so the iterations have converged. The simulator provides a timestep report which starts “Rep” meaning a report has been reached, and a “4” near the end of that line means that 4 non-linears were needed to converge.
The second timestep needs a total of 6 iterations. In this step we see that the aot at the last non-linear iteration drops from 9.00 to 0.08. This is a sign of “quadratic convergence” which will occur when the simulator is very close to the solution and reduces the aot by orders by magnitude at each non-linear iteration.

Details:

If you wish to tighten (reduce) the convergence criteria, there are minimum values below which SCONVP and SCONVS cannot be set. The pressure minimum is 0.01 and the effective saturation minimum 0.005. The pressure minimum can be ignored by setting the data item (1st in CVCRIT) negative. This behaviour is not documented, and may be changed in future releases.
The convergence tolerances are relaxed slightly at each non-linear iteration, to make convergence easier as the number of iterations increases. The hope is that this will allow a difficult timestep to be completed without significantly effecting the results. The factor used is:

Factor = 1.0 + “iteration number”/”maximum iterations”

An additional check is made that the sum of the residuals, which is a measure of total material balance error, is not excessive. This is a safety measure that can prevent convergence but rarely gets invoked in production cases.
Even if the solution change is too great and aot>1, the maximum residual may be small enough for the timestep to be converged. The criterion used is the variable SNLRMX set using the 11th data item in CVCRIT. The residual considered is the residual calculated using the previous iteration and is not the new residual resulting from the solution change. Because of this the test tends not to be effective very often.

Gain Option in ECLIPSE 300

“Gain” is an option, not used by default, which can significantly improve performance. The idea is to speed the code up by not taking an extra non-linear iteration if that iteration is likely to generate a very small solution change. We predict the behaviour of the next iteration based on the history of previous timesteps. We calculate a gain factor from previous timesteps and rather than using the aot we use:

EAOT=MIN(2.0*AOT*EGAIN,AOT)
Where: 
EAOT - Effective aot
AOT - Is the aot as calculated above
EGAIN - A measure of the expected improvement at the next step

The gain value is printed out with the non-linear debug.

The gain option can be switched on using the 68th switch of the OPTIONS3 keyword. The gain option can be controlled using the 19th switch of the OPTIONS3 keyword. If set to 1, the gain will not be used if any cell has changed state during the iteration. In this case you might expect the previous history to be less valid.
The option should be used with care. It can lead to dramatic changes in iteration sequence.

Tracking the source of the non-convergence problem

A single cell can cause non-convergence. As we increase the number of cells in a simulation, we increase the odds that a cell will cause non-convergence.

If the are only one or two cells in the reservoir model that are causing problems, we can identify these cells and check if there is any engineering or data reason which could explain why they are causing problems. For example the cells may be at or near well completions, in which case well control could be modified, or could be cells with very small pore volumes in which case the MINPV keyword could be used.

ECLIPSE 100 convergence

In ECLIPSE 100, you will get non-linear debug if you set “NEWTON=2” in RPTSCHED. This will produce output of the form:

IT= 0 CNV CELL MAT BAL DPRESS DSWAT DSGAS
OIL 1.00424( 28, 45, 3) 5.3D-03 0.00 0.00000 0.00000
WAT-0.00288( 9, 3, 3) -1.3D-07 0.00 0.00000 0.00000
GAS********( 5, 45, 1) -1.3D-02 0.00 0.00000 0.00000
LINIT= 5 NSCHP= 6 NCHOP= 0 NSTAT1,2,3= 50 5400 0 NTRAN= 321
IT= 1 CNV CELL MAT BAL DPRESS DSWAT DSGAS
OIL-1.99144( 5, 45, 1) 4.3D-02 -24.89 0.00026 -0.20000
WAT-0.16316( 2, 4, 4) -3.7D-06 -14.02 0.00490 0.00000
GAS********( 5, 45, 1) -3.1D-02 -24.89 0.00026 -0.20000
LINIT= 3 NSCHP= 195 NCHOP= 0 NSTAT1,2,3= 50 5370 30 NTRAN= 30
IT= 2 CNV CELL MAT BAL DPRESS DSWAT DSGAS
OIL-0.62319( 5, 45, 1) 1.7D-02 -21.15 0.00081 -0.01843
WAT-0.04162( 28, 5, 3) -1.3D-04 -30.47 0.00139 0.04000
GAS********( 5, 45, 1) -2.2D-02 -21.15 0.00081 -0.01843
LINIT= 3 NSCHP= 44 NCHOP= 3 NSTAT1,2,3= 50 5367 33 NTRAN= 25
IT= 3 CNV CELL MAT BAL DPRESS DSWAT DSGAS
OIL-0.30993( 5, 45, 1) -8.9D-05 -26.44 0.00088 -0.21134
WAT-0.04591( 28, 4, 3) -4.0D-05 -19.80 0.00666 0.11687
GAS********( 5, 45, 1) -1.4D-02 -26.44 0.00088 -0.21134

This shows the first 4 non-linear iterations (IT=0, IT=1, IT=2, IT=3) in a case that has convergence problems for the gas phase.
The first line shows IT=0, the first iteration, and CNV etc are column headers for the next 3 lines. The columns are:

CNV The worst residual for the OIL, WATer and GAS phases
CELL The cell that has the worst residual
MAT BAL The material balance for that cell, a measure of mass accuracy
DPRESS The change in pressure in that cell since the last iteration
DSWAT The change in water saturation since the last iteration
DSGAS The change in gas saturation since the last iteration

The residual for gas in all 4 iterations is shown as ******* which means that it is greater than the maximum printable value. It has a very high residual at each iteration for cell (5,45,1), so that is the cell that is causing problems.

After each iteration report above, the line starting LINIT= provides more information on what is happening within the model.

LINIT number of iterations required to solve the linearized equations.
NSCHP number of saturation changes that were altered to suppress possible oscillations.
NCHOP number of times the changes in P, Rs, or Rv were reduced to increase stability.
NSTAT1,2,3 is the number of cells in solution state 1,2,3 Solution state 1 means no oil is present in the cell Solution state 2 means both oil and gas are present in the cell Solution state 3 means no gas is present in the cell
NTRAN is the number of state transitions since the last non-linear iteration.

Any non-zero value of NSCHP or NCHOP increases material balance errors for the subsequent non-linear iteration and therefore reduces the chances of convergence.

Saturation chops

Some saturation chops can be avoided by adjusting relative permeability curves in such a way that the critical saturation is not the same as the lowest saturation value in the table. For instance, instead of

SWFN
0.2 0 7
0.3 0.07 4
0.4 0.15 3
0.5 0.24 2.5
0.6 0.33 2
0.8 0.65 1
0.9 0.83 0.5
1 1 0 /

try using:

SWFN
0.2 0 7
0.21 0 1*
0.3 0.07 4
0.4 0.15 3
0.5 0.24 2.5
0.6 0.33 2
0.8 0.65 1
0.9 0.83 0.5
1 1 0 /

The new saturation value at 0.21 may help convergence. It will not affect the initial fluids-in-place but will unfortunately slightly reduce the water mobility for water saturations between 0.2 and 0.3. This may be not so important to engineering accuracy.

Look for oscillations in the CNV for a phase. If one iteration has a positive value, the next iteration has a negative value, then the next is positive, then negative, etc. then there is perhaps a non-linearity in the system. These are sometimes associated with sudden changes in the slope of the relative permeability curves. If you have access to the SCAL program, you can plot these slopes and look for discontinuities. If you have access to a spreadsheet program then you can numerically calculate and plot the slopes. Remember that ECLIPSE will use all the values of saturation and relative permeability that you give in the table without any smoothing. You should therefore try to avoid tables such as:

SWFN
0.2 0 7
0.21 0 1*
0.3 0.07 4
0.301 0.07 4
0.398 0.14 3
0.4 0.15 3
0.401 0.17 3
0.402 0.19 3
0.5 0.24 2.5
0.6 0.33 2
0.8 0.65 1
0.9 0.83 0.5
1 1 0 /

The table above has saturation values that are too close to each other and the slopes of the relative permeabilities shows severe changes. You should also try to avoid tables such as:

SWFN
0.2 0 7
0.3 0. 4
0.5 0.01 2.5
0.51 0.60 2
0.8 0.68 1
0.9 0.83 0.5
1 1 0 /

The table above has a very sudden krw change from krw=0.01 at Sw=0.5 to krw=0.60 at Sw=0.51 and will certainly cause convergence problems.

TUNINGDP output

The NEWTON switch in RPTSCHED will produce extra information in the case of TUNINGDP:

PTRG is the target pressure change; the default is 1 psi in Field units.
STRG is the target saturation change; the default is 0 if TUNINGDP is not used and 0.01 if it is used.
MDDP is the maximum Pressure change for convergence.
MDDS is the maximum Saturation change for convergence.

If you use TUNINGDP:

ECLIPSE will solve the linear equations to a tighter tolerance
the convergence is reached if either the residual or the solution change criteria is small enough.

If you don‟t use TUNINGDP, only the residual is used to test for convergence.

Output of reason for non-linear failure

If you set DEBUG item 1 in ECLIPSE 100 to be > 1, ECLIPSE will output a number that shows the reason that a non-linear iteration has failed to converge. The 5 possible values are: 1 - Exceeds maximum tolerable pressure change DDPLIM
2 - Exceeds maximum tolerable saturation change DDSLIM
3 - Exceeds maximum tolerable material balance RSUM
4 - Exceeds tolerable for maximum norm RNMAX
5 – Uses but fails the TRGDDP/TRGDDS test
DDPLIM and DDSLIM are described in record 3 of the TUNING keyword.

A value of 3 means that the normalised residual is greater than the maximum allowed. That allowed maximum is a linear combination of items 3 and 7 of record 2 of the tuning parameter. The maximum residual also depends on which newton iteration we are on.

A value of 4 means that the normalised residual is greater than the maximum allowed. That allowed maximum is a linear combination of items 2 and 6 of record 2 of the tuning parameter. The maximum residual also depends on which newton iteration we are on.

A value of 5 means that the change in the solution is less than TRGDDP and TRGDDS. These are the minimum values of DDPLIM and DDSLIM until that point in the simulation. This is used in cases where there is high throughput.

ECLIPSE 300 convergence

In ECLIPSE 300, you can have either a visual display of the grid blocks causing convergence problems or you can look at numerical output in the same way as for ECLIPSE 100. Starting with the 2009.1 version, ECLIPSE 100 may also have the visual display option.

Adding CONV to the RPTRST keyword will send two new outputs to the restart files. Each cell will have two new variables that will be used to count the number of times that cell has been one of the most difficult cells to converge. At the beginning of the simulation each cell will have its counter set to zero. At every timestep, the 10 most difficult cells will have their counter increased by 1. At the end of the run you can display the cells with the most problems.

Visual option in RPTRST – ask for CONV

Output of cells which are causing convergence problems. By default, CONV=10 is set so that the worst 10 cells will be output

Output:

CONV_VBR: Worst cells based on volume balanced residual
CONV_PRU: Worst cells based on pressure updates

Alternatively, you will get non-linear debug if the 8th data item in the DEBUG3 keyword is set greater than 0. Typical output is of the form:

Iteration 0 linears req 7
DX Pressure 0 -40.075969 25 32 4 F 1.469590
DX Comp 1 -0.000375 25 32 1 T 0.010000
DX Comp 2 -0.000980 25 32 1 T 0.010000
DX Comp 3 -0.025419 25 32 1 F 0.010000
DX Comp 4 -0.002318 25 32 1 T 0.010000
DX Comp 5 -0.000318 25 32 3 T 0.010000
DX Comp 6 0.009711 25 32 1 T 0.010000
DX Comp 7 0.008562 25 32 1 T 0.010000
DX Comp 8 0.004396 25 32 1 T 0.010000
DX Comp 9 0.001465 25 32 1 T 0.010000
DX Comp 10 0.000907 25 32 4 T 0.010000
NLStep= 0 lin= 7 aot= 27.27 Rmax=0.8514E+00 Rsum=0.2739E-03 egain=0.2624E-01
Iteration 1 linears req 5
DX Pressure 0 -0.563835 25 32 1 T 1.592056
DX Comp 1 0.000035 25 32 1 T 0.010833
DX Comp 2 0.000063 25 32 1 T 0.010833
DX Comp 3 0.001887 25 32 1 T 0.010833
DX Comp 4 0.000242 25 32 1 T 0.010833
DX Comp 5 0.000111 25 32 1 T 0.010833
DX Comp 6 0.000270 26 32 1 T 0.010833
DX Comp 7 0.000238 26 32 1 T 0.010833
DX Comp 8 0.000127 26 32 1 T 0.010833
DX Comp 9 0.000044 26 32 1 T 0.010833
DX Comp 10 -0.000317 25 32 2 T 0.010833
NLStep= 1 lin= 5 aot= 0.38 Rmax=0.1921E-01 Rsum=0.1448E-04 egain=0.3165E-01
Max changes:pres 40.6 25 32 4 temp 0.00 0 0 0
oil satn 0.516E-01 17 7 1 gas satn -0.178E-01 17 8 1
wat satn -0.939E-03 25 32 4 eng dens 0.00 0 0 0
Throughput ratio:avrg 0.404E-01 max 0.192 26 32 2
MIF ; 103.0 9.00 6.9094 .01725 8973.3 157.46 62000. 3531.4 0.0 60500. 2 2%

This output shows two non-linear iterations leading to a timestep report.

The first line:

Iteration 0 linears req 7

tells us that the iteration 0 (the initial estimate) needed 7 linear iterations to solve the linear problem.

The next 11 lines consist of one line for each of the solution variables showing the largest change in that variable in any cell during this iteration. The solution variables for each grid block are the pressure in the grid block and the molar density of each hydrocarbon component, and a water term. This model has 9 hydrocarbon components. Water is written as component 10.

DX Pressure 0 -40.075969 25 32 4 F 1.469590
DX Comp 1 -0.000375 25 32 1 T 0.010000
DX Comp 2 -0.000980 25 32 1 T 0.010000
DX Comp 3 -0.025419 25 32 1 F 0.010000
DX Comp 4 -0.002318 25 32 1 T 0.010000
DX Comp 5 -0.000318 25 32 3 T 0.010000
DX Comp 6 0.009711 25 32 1 T 0.010000
DX Comp 7 0.008562 25 32 1 T 0.010000
DX Comp 8 0.004396 25 32 1 T 0.010000
DX Comp 9 0.001465 25 32 1 T 0.010000
DX Comp 10 0.000907 25 32 4 T 0.010000

In each of these lines, DX means the solution change. The first line is the Pressure change. The largest pressure change was an increase of 40.075969 psi in cell (25,32,4), which happen to contain an injecting completion. The “F” on that line means “False”, in that the pressure variable has not converged, as the pressure change is greater than 1.469590, which is the maximum pressure change, allowed for convergence for this iteration.
The second DX line shows the largest change in the molar density, expressed as a saturation equivalent, for component 1. This increase of 0.000375 was in cell (25,32,1) and is less than the convergence maximum of 0.01, so that the component 1 variable is considered to be converged. In fact all the components have converged except for component 3.
The non-linear iteration however has not converged since two of the variables (pressure and component 3) are not yet converged.

The next line is a summary of the first iteration (iteration 0).

The line:

Iteration 1 linears req 5

is the start of the report on the second non-linear iteration, Iteration 1, which needed 5 linear iterations.

The next 11 lines are the solution changes. They are similar to the reported changes for the first non-linear iteration except that now both the pressure and component 3 have changed by less than the new convergence criteria. The non-linear iterations have now converged.

NLStep= 1 lin= 5 aot= 0.38 Rmax=0.1921E-01 Rsum=0.1448E-04 egain=0.3165E-01

is the report of the converged iteration showing that the maximum residual at the beginning of that iteration was down to 0.1921E-01 and the sum of residuals was down to Rsum=0.1448E-04.

There then follows a report on the changes during that timestep:

Max changes:pres 40.6 25 32 4 temp 0.00 0 0 0
oil satn 0.516E-01 17 7 1 gas satn -0.178E-01 17 8 1
wat satn -0.939E-03 25 32 4 eng dens 0.00 0 0 0

The maximum pressure and saturation changes are reported, as well as the cells in which this change occurred. In the case of thermal runs, the maximum temperature and energy density changes are also reported.
The maximum throughput is reported next. Throughput is defined as the volume flowing through a cell divided by the pore volume of the cell. If the throughput is too high (say higher than 0.5) it could cause convergence problems, and the pore volume of the cell that has the high throughput should be examined.

 Throughput ratio:avrg 0.404E-01 max 0.192 26 32 2

The last of these lines is the report of production, etc. for the timestep:

 MIF ; 103.0 9.00 6.9094 .01725 8973.3 157.46 62000. 3531.4 0.0 60500. 2 2%

Moving goalposts

ECLIPSE 300 non-linear convergence criteria have 'moving goalposts'. The convergence tolerances are relaxed slightly at each non-linear iteration. Effectively, convergence becomes easier as the number of iterations increases.

The idea is that if you have reached the maximum number of iterations (call that Nmax) and you are close (within a factor of two) to the convergence criteria, then you don't want to chop the timestep and waste all the work you have done so far. So the criteria relaxes by 1/Nmax every Newton iteration.

An example is shown below. The maximum number of non-linear iterations is 12, the units are metric, and the first Newton iteration uses the default convergence criteria (0.1 atm pressure and 0.01 for component specific volume). For later Newton iterations, these criteria were relaxed by 8.33% (=1/12) with each Newton iteration. After 12 Newton iterations, the criteria are doubled to 0.2 atm pressure and 0.02 for component specific volume.

Iteration 0 linears req 4 NTOTNL 4076 
DX Pressure 0 2.803147 38 16 10 Ind:GLOB F 1.469590 
DX Comp 1 0.000262 34 31 6 Ind:GLOB T 0.010000 
DX Comp 2 0.000042 34 31 6 Ind:GLOB T 0.010000 
.. 
NLStep= 0 lin= 4 aot= 1.91 Rmax=0.1149E-02 Rsum=0.8487E-05 egain=0.9926E+00 
DCHOP2: 1 cells chopped, Try= 1 
Iteration 1 linears req 7 NTOTNL 4077 
DX Pressure 0 -2.799664 38 16 10 Ind:GLOB F 1.592056 
DX Comp 1 -0.003126 38 16 10 Ind:GLOB T 0.010833 
DX Comp 2 -0.000239 38 16 10 Ind:GLOB T 0.010833 
.. 
NLStep= 1 lin= 7 aot= 1.91 Rmax=0.3573E-01 Rsum=0.1269E-05 egain=0.3326E+01 
DCHOP2: 1 cells chopped, Try= 1 
Iteration 2 linears req 4 NTOTNL 4078 
DX Pressure 0 2.694699 38 16 10 Ind:GLOB F 1.714522 
DX Comp 1 -0.048358 38 16 10 Ind:GLOB F 0.011667 
DX Comp 2 -0.002969 38 16 10 Ind:GLOB T 0.011667

Discussion on Non-Linearity

A few general problems are:

In explicit cells, the algorithm can try to extract more fluid from a grid block than is present in the grid block. This happens because flows are calculated from the solution of the previous timestep, and not from the current timestep. The only thing that can be done is to reduce the timestep. This will only be an issue if you using the IMPES or AIM solution methods. The default in ECLIPSE 100 is FULLIMP, but the default in ECLIPSE 300 is AIM unless you are using some options such as Radial, Dual Porosity or Thermal.
Flow reversals are a major non-linearity.
The group control algorithm sometimes wants to change well rates at every non-linear iteration to reflect the latest calculated conditions in the reservoir. Changing well rates every non-linear iteration can however lead to poor convergence. The simulators therefore by default only recalculate group control parameters for the first 4 non-linear iterations, then keep the group controls unchanged for the remaining non-linears. This value of 4 can be changed using NUPCOL, but we recommend leaving the value unchanged.
Non-monotonic VFP tables can cause convergence problems. VFP tables are always checked for monotonicity in ECLIPSE 300, and the check can be switched on in ECLIPSE 100 by using the EXTRAPMS keyword.

Non-Linear Divergence

Sometimes oscillation or “flip-flop” occurs. This is can be seen in the debug output:

LINIT=10 NSCHP= 6 NCHOP= 0 NSTAT1,2,3= 0 1793 1213 NTRAN= 3
IT= 6 CNV CELL MAT BAL DPRESS DSWAT DSGAS
OIL -0.03205( 10, 11, 3) 1.1D-06 -84.80 0.00002 -0.01805
WAT-0.00004( 10, 11, 3) -2.6D-08 -84.80 0.00002 -0.01805
GAS-0.07982( 13, 10, 3) -3.5D-06 -82.33 0.00006 -0.12626
LINIT= 6 NSCHP= 0 NCHOP= 0 NSTAT1,2,3= 0 1793 1213 NTRAN= 0
IT= 7 CNV CELL MAT BAL DPRESS DSWAT DSGAS
OIL -0.00004( 10, 18, 3) 5.0D-09 -84.65 0.00006 -0.00989
WAT 0.00000( 8, 17, 3) -1.0D-13 -86.77 0.00006 -0.04626
GAS 0.00354( 13, 10, 3) 2.0D-09 -81.34 0.00006 -0.11912
LINIT=10 NSCHP= 6 NCHOP= 0 NSTAT1,2,3= 0 1793 1213 NTRAN= 3
IT= 8 CNV CELL MAT BAL DPRESS DSWAT DSGAS
OIL -0.03205( 10, 11, 3) 1.1D-06 -84.80 0.00002 -0.01805
WAT-0.00004( 10, 11, 3) -2.6D-08 -84.80 0.00002 -0.01805
GAS-0.07982( 13, 10, 3) -3.5D-06 -82.33 0.00006 -0.12626
LINIT= 6 NSCHP= 0 NCHOP= 0 NSTAT1,2,3= 0 1793 1213 NTRAN= 0
IT= 9 CNV CELL MAT BAL DPRESS DSWAT DSGAS
OIL -0.00004( 10, 18, 3) 5.0D-09 -84.65 0.00006 -0.00989
WAT 0.00000( 8, 17, 3) -1.0D-13 -86.77 0.00006 -0.04626
GAS 0.00354( 13, 10, 3) 2.0D-09 -81.34 0.00006 -0.11912

The solution for non-linear iteration 6 is the same as that for non-linear iteration 8, and the solution for non-linear iteration 7 is the same as that for non-linear iteration 9.

Non-Linear Equation Convergence Failure - Case Study

The problem when I input Pc=0 in the model amd plot saturation versus distnace the is not any shock fron. Then when I input the Pc with value for the model and plot sw versus distance It gave shock front for the firs step. the step for the model 600*0.005.

The message looks like this:

@--PROBLEM AT TIME 2.3 HOURS ( 1-JUL-2015):
@ NON-LINEAR EQUATION CONVERGENCE FAILURE 
@ ITERATION LIMIT REACHED - BUT TIME STEP 
@ ACCEPTED BECAUSE IT IS TOO SMALL TO CHOP

Why?

The simulator is reporting to be unable to solve the material balance equations within the accuracy desired by the user and within the max. number of Newton (thus nonlinear, outer) iterations, even when it decreases the simulation timestep to the minimum value, so it accepts the solution at that minimum timestep length and proceeds further (see the TUNING keyword; for dual porosity/permeability and other 'high throughput" runs see also TUNINGDP). The error of this type is generally 'unacceptable' (as opposite to not fully converged solutions of the ('inner loop') linear systems at every Newton iteration, which is not necessarily a problem as long as the material balance, i.e., Newton iteration is solved within the desired accuracy; non-convergent linear iterations may lead to slower runs and issues with nonlinear iterations later in the run, however).

This type of convergence error is generally throughput related, put simply it means that the volume of fluid passing from one grid cell to the grid cell adjacent to it is to large. This can be caused by varying pore volumes but more likely due to pinched out grid cells which have peculiar angles as a rsult of faulting implemented into the model. You can create properties to make the pinched out cells inactive which will help avoid this throughput convergence.

This kind of convergence failure is related to INPUT data and/or modeling:
1) EXTRAPOLATION in PVT or VFP table. Use the EXTRAPMS keyword to check.
2) Initialization: 2 different equilibrium region without sealing fault.
...
There are some other causes related to DUALPORO/

Possible solutions

If the simualtor failed to solve the equation within the allocated timestep - to a degree this is accceptable, it will affect the accuracy of the results.
If model is full of convergence errors like this you really need to analyse whats causing them.

If the model run within acceptable CPU time you can ignore the warning. Almost all the models have this kind of warning. It's not necessary to avoid it.
However, if your mode run very slow and have convergency problems, don't try to tuning the "default parameters of the solver", which is reasonable value defined by the developer. Otherwise, focus on the model itself, including the geological properties,PVT, rock and the production data. For example, reservoir model is a person and simulator is shoes. if the person run slow, the problem probably from the physic body, nor the shoes. any method to fix the shoes may not work.

Check if model has too many LGRs or too many cells where the model pinches out.

Check out there are certain keywords for non linear equations.

Use EXTRAPMS to see if there is any extrapolation in PVT table.

Correct the well PI by using WPIMULT keyword.

Change time stepping control by using the TUNING keyword or by assigning very small report steps using TSTEP

Load your SCAL data in SCAL software and check if there is very sharp change in curves derivative

A simple workaround is to enter the NSTACK keyword into the Runspec section of your ECLIPSE data set, followed by a number like 50 or 100, this will increase the amount of cpu memorary allocated to solving the equation just write it

NSTACK 
50/

One more thing can be done is increase the number of linear iterations that ECLIPSE will use to solve the equation, this gives it a better chance to get the correct answer and will therefore maintain the accuracy of your results. Beware this can either cause simulation times to speed or or slow down depending on the model. You should not set NSTACK higher than whatever you set LITMAX to. To set the LITMAX (linear iterations limit) in the SCHEDULE section write

TUNING 
/ 
/ 
2* 50/

Be aware of setting the limiting the NSTACK and LIYMAX according to the cells sizes and the number, these if over used may drastically change your result.

To see cells causing convergence problems RPTRST use keyword with the Mnemonic CONV, this will allow the output of cells which are causing convergence problems.
Setting CONV=n focusses on the n worst cells. By default, n=1(E100) and n=10(E300). Different things has mentioned above can cause convergence problems.

Try MINPV with PINCH Keyword to limit convergence issues.

Try and increase LITMAX. You can try also to increase the max number of non linear iterations(check 3rd record of the Tuning keyword).

Deactivate unnecessary cells like shale.

- generally to avoid numerical difficulties the Grid must be regular and orthogonal,
- eliminate cells with low pore volume(MINPV),
- avoid too large ratio of Pore Vol between 2 neigbour cells,
- avoid cells with low PV and large T.
- LGR are time consuming use locally, do not locate wells close to LGR boundary.
- record 3 of TUNING keyword: Decrease item1, increase item 3 increase item5 and item 6.
- PVT and SCAL data must be smooth. Plot and confirm same.
- WVFPEXP: will help if convergence is due to well control in THP.
- confirm there is no communication between different PVT, Equilibration regions.
- after initialisation, check the stability with the well shut in Run for a year+. Plot in the SUMMARY section MAXDPR, MAXDSW, MAXDSG.

You can follow the suggestions above but sometimes is is anavoidable.
In that case use in RUNSPEC:

MESSAGES 
9* 100000/

Source

The theme was discussed here http://lnkd.in/KEuFRc

Convergence problem Checklist

How to identify cells with convergence problems

Add CONV to RPTRST
To get more information for ECLIPSE 100, add „NEWTON=nn‟ to RPTSCHED
To get more information for ECLIPSE 300, add:

DEBUG3
 7* 1 /

What to look for in different Sections

Basic
Check all problem and warning messages. See Eclipse Severity Levels

GRID data

Check whether you need to set a minimum pore volume.

MINPV

Use MINPV to remove cells with small pore volumes. If a full-field model has convergence problems and does NOT have a minimum pore volume, adding MINPV is the most likely change that will improve performance.
Take care to choose a suitable value of MINPV that does not significantly change the total pore volumes in each region, and that does not remove high-permeability streaks or thin shale layers within the reservoirs. If MINPV is being used to remove pinch-outs, also use PINCH to connect across those thin cells that represent pinch-outs. Within these limitations, using MINPV should both improve performance and give unchanged production results.

LIMITS

Look for the table headed LIMITS in the ECLIPSE 300 PRT file.
This is a table of minimum and maximum values for each GRID array. Check the ratio of Maximum/Minimum values at the left of output, and investigate if any of these ratios is very large.

Grid Array Maximum value Minimum value Non-zero minimum Ratio max/min
---------- ---------------------------- ---------------------------- ---------------------------- --------------
TOPS 7571.0 ( 1 5 9) 7021.0 ( 1 1 1) 7021.0 ( 1 1 1) 1.078336
MIDS 7926.0 ( 1 5 9) 7146.0 ( 1 1 1) 7146.0 ( 1 1 1) 1.109152
DEPTH 7926.0 ( 1 5 9) 7146.0 ( 1 1 1) 7146.0 ( 1 1 1) 1.109152
DX 220.0 ( 1 1 1) 55.0 ( 1 6 2) 55.0 ( 1 6 2) 4.0
DY 220.0 ( 1 3 1) 55.0 ( 1 1 1) 55.0 ( 1 1 1) 4.0
DZ 710.0 ( 1 1 3) 250.0 ( 1 1 1) 250.0 ( 1 1 1) 2.84
DZNET 710.0 ( 1 1 3) 250.0 ( 1 1 1) 250.0 ( 1 1 1) 2.84
NTG 1.0 ( 1 1 1) 1.0 ( 1 1 1) 1.0 ( 1 1 1) 1.0
PORO 0.10 ( 1 1 1) 0.15E-01 ( 1 3 4) 0.15E-01( 1 3 4) 6.666667
PORV 459036.7 ( 1 1 3) 323.26 ( 1 1 4) 323.26 ( 1 1 4) 1420.0
PERMX 300.0 ( 1 1 4) 0.0 ( 1 1 1) 1.5 ( 1 1 6) 200.0
PERMY 300.0 ( 1 1 4) 0.0 ( 1 1 1) 1.5 ( 1 1 6) 200.0
PERMZ 7.5 ( 1 1 4) 0.0 ( 1 1 1) 0.75 ( 1 1 6) 10.00
TRANX 8.45 ( 1 1 4) 0.0 ( 1 1 1) 1.20383 ( 1 1 6) 7.042254
TRANY 8.45 ( 1 1 4) 0.0 ( 1 1 1) 1.20383 ( 1 1 6) 7.042254
TRANZ 1.49 ( 1 1 4) 0.0 ( 1 1 1) 0.1105796( 1 1 5) 13.45455
MULTX 1.0 ( 1 1 1) 1.0 ( 1 1 1) 1.0 ( 1 1 1) 1.0
MULTY 1.0 ( 1 1 1) 1.0 ( 1 1 1) 1.0 ( 1 1 1) 1.0
MULTZ 1.0 ( 1 1 1) 0.0 ( 1 1 9) 1.0 ( 1 1 1) 1.0
ROCKV 6028682. ( 1 1 6) 1939592. ( 1 4 4) 1939592. ( 1 4 4) 3.108222
MULTPV 1.0 ( 1 1 1) 1.0 ( 1 1 1) 1.0 ( 1 1 1) 1.0
MINPVV 0.1E-05 ( 1 1 1) 0.1E-05 ( 1 1 1) 0.1E-05 ( 2 1 1) 1.0
------------------------------------------------------------------------------------------------------------------------------

In-Place LGRs

By default, ECLIPSE 100 will solve LGRs using local time stepping.
With local time stepping, the global grid can take a large time step that is not limited by the size of the time step needed by the LGR.

The local time stepping method is semi-explicit and therefore is potentially unstable. The alternative is to use the LGRLOCK to “lock” the timesteps together and turn off the local time stepping, so that the global timesteps are the same as the local timesteps. With LGRLOCK, the method is fully implicit and unconditionally stable.
The keyword LGRFREE will allow local timestepping to start again, if the convergence issue only occurs for a few timesteps.

SCAL data

Check the slopes of all relative permeability tables.

Example:

-- SOIL KROW KROG
SOF3
0.181 0 0
0.283 0.0001 0.0001
0.385 0.0015 0.0015
0.436 0.0124 0.0124
0.483 0.0217 0.0217
0.588 0.0939 0.0939
0.686 0.3499 0.3499
0.689 0.3501 0.3501
0.761 0.7323 0.7323
0.837 0.9887 0.9887
0.863 0.9978 0.9978
0.879 1 1
0.880 1 1 /

The curve is smooth and looks reasonable.
The Soil and krow column see to be correct. The Soil column is monotonically increasing, and has values between 0 and 1. The krow column is also monotonically increasing, and has values between 0 and 1.

Try to plot the derivative d(krow)/d(Soil). It is expected the derivatives to be smoothly varying.

But in fact, there are sharp changes in slope around Sw = 24% and around Sw = 31%. These value correspond to Soil = about 68% and 76%. At around 68% oil there are two oil saturations with very close values, 0.686 and 0.689, that have very similar krow values of 0.3499 and 0.3501.

The SOF3 table does not need saturation values that close together, and in fact very closely spaced values in any table may cause ECLIPSE to do extra computing work as it has to find exactly which two point to interpolate between. In this case we can fix the problem by removing either one of the two entries.

Removing very small kr values

We can however also make an engineering improvement to the table. Consider the entry at Soil=0.283.

-- SOIL KROW KROG
SOF3
0.181 0 0
0.283 0.0001 0.0001
0.385 0.0015 0.0015
0.436 0.0124 0.0124
0.483 0.0217 0.0217
0.588 0.0939 0.0939
0.686 0.3499 0.3499
0.689 0.3501 0.3501
0.761 0.7323 0.7323
0.837 0.9887 0.9887
0.863 0.9978 0.9978
0.879 1 1
0.880 1 1 /

That line in the table says that at 28.3% oil saturation the relative permeability of oil to water and to gas are both 1E-4. It is possible but unlikely that this 1E-4 is a value that has been measured experimentally. From a simulation point of view it makes more sense to set that value to 0. ECLIPSE has special code (known as the Appleyard Chop) that speeds up convergence at the point when a phase first becomes mobile. That code is wasted if the change in relative permeability is from 0 to 1E-4, but is very effective if the change is from 0 to a higher value.

Critical saturations

Consider a grid block with pore volume of 100 m3. Assume it contains 80 m3 of gas and 20 m3 of water. The water saturation is 20%.

Suppose:

The pore volume decreases by 1% due to reservoir pressure changes
There is no flow in or out of that grid block.

The gas will compress but water is relatively incompressible, so we will have approximately 79 m3 of gas and 20 m3 of water. The water saturation will increase from 20/100 = 20% to 20/99 = 20.2%

As you produce/inject, the pressure in every block will change.
Therefore the pore volume of each grid block will change.
If the initial water saturation (SWL) in each grid block was 20%, Sw will now be >20% in many blocks.
If SWCR=SWL=20% then all those blocks have mobile water even if Sw is only 20.0001%. This water WILL flow.
Setting SWCR=SWL+0.01% will stabilize the model without changing the OOIP.

The same applies to the critical saturation of any phase. The original SOF3 table is

-- SOIL KROW KROG
SOF3
0.181 0 0
0.283 0.0001 0.0001
0.385 0.0015 0.0015
0.436 0.0124 0.0124
0.483 0.0217 0.0217
0.588 0.0939 0.0939
0.689 0.3501 0.3501
0.761 0.7323 0.7323
0.837 0.9887 0.9887
0.863 0.9978 0.9978
0.879 1 1
0.880 1 1 /

As discussed above, we could stabilize the oil by changing the 1E-4 values to zero. If for some reason we do not wish to do that, we could instead add a relative permeability of 0 at an oil saturation of 18.2%:

-- SOIL KROW KROG
SOF3
0.181 0 0
0.182 0 0
0.283 0.0001 0.0001
0.385 0.0015 0.0015
0.436 0.0124 0.0124
0.483 0.0217 0.0217
0.588 0.0939 0.0939
0.689 0.3501 0.3501
0.761 0.7323 0.7323
0.837 0.9887 0.9887
0.863 0.9978 0.9978
0.879 1 1
0.880 1 1 /

SCAL DATA: other issues

Hysteresis

Does the model run better without hysteresis?
Are the drainage and imbibition end-points consistent?

End-point scaling

Are the end-point values extreme (for instance SWL=>1)?
Are the end-points too close?
If you are using SCALECRS, is there a severe slope discontinuity at both ends of the 2-phase mobile regions? For example in an oil-water problem, do you have a discontinuity when (1-Sowcr) approaches SWU? In this case, try using SCALELIM to limit the value of (1-Sowcr) used in scaling.

PVT data

Total Compressibility Check

An example:

@--WARNING AT TIME 0.0 DAYS (19-OCT-1982):
@ NEGATIVE COMPRESSIBILITY FOUND IN GAS
@ PRESSURE TABLE 1 AND OIL PRESSURE
@ TABLE 1 AT A SAMPLE PRESSURE VALUE
@ 290.56207 . ADJUST SATURATED FLUID
@ PROPERTY VALUES AT THIS PRESSURE.
@ NEGATIVE COMPRESSIBILITIES OCCUR
@ FOR GAS SATURATIONS LESS THAN 1.00000

Although this is only a warning, it is probably NOT safe to ignore the warning.

In black oil models, both ECLIPSE 100 and ECLIPSE 300 check for positive compressibility of each single reservoir fluid as the PVT data is read (the formation volume factor must be a monotonically decreasing function of pressure assuming all other variables are held fixed).
ECLIPSE also checks that mixtures of saturated oil and gas have a positive total compressibility even when there is mass transfer between the two phases. For example, increasing the pressure of a cell that contains both oil and gas will:

transfer some gas from the gas phase to the oil phase
swell the oil due to extra dissolved gas
decrease the remaining gas volume due to compression

The oil volume will therefore increase with increasing pressure, and the gas volume will decrease with increasing pressure. The total (oil+gas) volume must however decrease. This will only happen if the reduction in the volume of gas is greater than the increase in the volume of oil.

The total compressibility of the oil and gas system is given by the expression:

Ct = Sg*Ctg + So*Cto 
Ctg = [-dBg/dP + dRv/dP * (Bo-Rs*Bg) / (1-Rs*Rv)]/Bg 
Cto = [-dBo/dP + dRs/dP * (Bg-Rv*Bo) / (1-Rs*Rv)]/Bo

Where Sg, So, Bg,Bo, Rs, and Rv have their usual meaning, and the pressure derivatives are taken along the saturation curve.

In the special case of Rv=0, this simplifies to: 
Ctg = [-dBg/dP ] / Bg 
Cto = [-dBo/dP + dRs/dP * Bg] / Bo

For each PVT region in the simulation grid a pressure range is selected from the corresponding oil and gas PVT data tables to span the complete range of pressure data in the two tables. This pressure range is then subdivided into 30 equally spaced pressure nodes to evaluate the total hydrocarbon compressibility.

At each pressure node two limiting compressibility values are calculated corresponding to Sg=0 and Sg=1:

Ct = Ct,o at Sg = 0, So = 1 
Ct = Ct,g at Sg = 1, So = 0

If either or both values are negative at a pressure node, you will get a warning message.

The highest pressure used is the maximum of the oil pressure (maximum bubble point pressure in PVTO or maximum pressure specified in PVDO) and the gas pressure (maximum dew point pressure specified in PVTG or maximum gas pressure specified in PVDG). The maximum sample pressure at which total compressibility is checked can be overridden by choosing the second item in the PMAX keyword.

If the range of sample pressures extends above the maximum bubble point specified in the PVTO table or above the maximum dew point specified in the PVTG table, then ECLIPSE is forced to extrapolate above the highest entered bubble point or dew point (see below). In this case, it is not unlikely that negative compressibilities could occur as a result of extrapolation.

In ECLIPSE 100 only, if the 21st switch in the DEBUG keyword is set > 0, ECLIPSE 100 reports the total gas and oil compressibilities and as well as the saturated values of Rs, Bo and Rv,Bg for the sample range of pressures to the debug file (extension .DBG) in tabular form. Whenever any negative compressibilities are encountered in the data checking, this debug output is automatically switched on. You can use this information to diagnose any negative compressibilities reported by ECLIPSE 100.

Typical output from DEBUG(21)>0 is:

Total compressibility report:
Tcoil = total compressibility with Sg=0
Tcgas = total compressibility with So=0
Tc = Sg.Tcgas + So.Tcoil
Minimum oil pressure 0.00000
Maximum oil pressure 8.23759
Minimum gas pressure 0.00000
Maximum gas pressure 196.1330
=============================================
Press Rssat Rvsat Bg Bo DBg/DPo DBo/DPo DRs/DPo DRv/DPo Tcoil Tcgas
10.0000 1.9000 0.0000 1.0150 1.0230 -1.0502 0.0254 0.2039422 0.000000 0.1774885 1.0346498
16.7632 2.1424 0.0000 0.1336 1.0395 -0.0152 0.0009 0.0254929 0.000000 0.0023769 0.1135191
23.5264 2.3148 0.0000 0.0755 1.0459 -0.0058 0.0009 0.0254929 0.000000 0.0009356 0.0772953
30.2896 2.4872 0.0000 0.0493 1.0523 -0.0029 0.0010 0.0254929 0.000000 0.0002841 0.0589455
37.0528 2.6596 0.0000 0.0352 1.0588 -0.0015 0.0010 0.0254929 0.000000 -0.0000670 0.0421443
43.8160 2.8320 0.0000 0.0272 1.0654 -0.0009 0.0010 0.0254929 0.000000 -0.0002702 0.0343437
50.5792 3.0044 0.0000 0.0213 1.0721 -0.0012 0.0010 0.0254929 0.000000 -0.0004215 0.0543732
57.3424 3.1768 0.0000 0.0155 1.0400 -0.0006 0.0009 0.0254929 0.000000 -0.0005183 0.0397541
64.1057 3.3493 0.0000 0.0142 1.0401 -0.0001 0.0009 0.0254929 0.000000 -0.0005506 0.0045428
70.8689 3.5217 0.0000 0.0138 1.0401 -0.0001 0.0009 0.0254929 0.000000 -0.0005610 0.0044073

Since the total compressibility equation is:

Ct = Sg*Ctg + So*Cto

Then when P=70.8689 we can substitute the values of Ctg=Tcoil anf Ctg=Tcgas:

Ct = Sg*0.0044073 – So*0.0005610

and Ct will be negative when So is close to 1 and Sg is close to 0.

We have in this case of Rv=0, so the equations for Ctg and Cto are:

Ctg = [-dBg/dP ] / Bg 
Cto = [-dBo/dP + dRs/dP * Bg] / Bo

To correct the data in this case we can either decrease the absolute value of dBo/dP or increase the value of DRs/dP or of Bg.

PVT DATA: Extrapolation

An example:

@--WARNING AT TIME 0.0 DAYS (19-OCT-1982):
@ NEW SOLUTION OBTAINED BY EXTRAPOLATION OF
@ BOTH OIL AND GAS PVT TABLES

This is only a warning, but it may NOT be safe to ignore it.

In some cases, if the pressure is extrapolated severely, then the linearly-extrapolated formation volume factors may become unphysical. It is recommended that the highest bubble point nodes in the PVTO table are constructed so avoid extrapolations above the highest Rs entered in the PVTO table during the simulation (similarly, in a run with vaporized oil present, it is recommended that the highest dew point nodes in the PVTG table are constructed to avoid extrapolations above the maximum Rv entered during the simulation).

The keyword EXTRAPMS asks ECLIPSE 100 to warn you whenever extrapolations are made to PVT (or VFP) tables. In ECLIPSE 100, ALWAYS add keyword:

EXTRAPMS
4 /

ECLIPSE 300 will warn you by default
If insufficient PVT data is supplied, ECLIPSE may extrapolate the PVT table data to inaccurate or non-physical values.

Question:
If your initial reservoir pressure is 4800 psia, is it OK to have PVT tables with a maximum pressure value of 4800 psia?
Answer:
No, because any injection is likely to be at a pressure greater than 4800 psia.

If you are checking extrapolated values yourself, you should know that extrapolation is not linear in FVF or Viscosity. ECLIPSE stores PVT tables internally as the reciprocals of FVF and Viscosity* FVF.

With EXTRAPMS value of 4, you may get extra information:

@--WARNING AT TIME 0.0 DAYS (19-OCT-1982):
@ NEW SOLUTION OBTAINED BY EXTRAPOLATION OF
@ BOTH OIL AND GAS PVT TABLES
Extrapolation in PVTO table number 1
Po-Pbub above highest value on undersaturated line
Rs = 1.27027 MSCF/STB Po = 4820.658 PSIA
Pbub = 4014.700 PSIA
Cell (20,20,7) State 3
Sw = 0.12000 Sg = 0.00000 So = 0.88000

SOLUTION data

Check that the initial solution is stable

If you are using EQUIL, it should be in equilibrium
For enumerated solutions, it depends on what you want to do
For RESTARTs, it won‟t be usually be in equilibrium

To check for stability:

Run the model with no wells and no aquifers
Start with small timesteps and increase the timestep size
Check that there are no convergence problems
Check that field pressures and fluids-in-place do not change Plot SUMMARY quantities such as FPR, FOIP, …

Verify that it doesn't have saturation or pressure variations. It's easy and should always be done before proceeding with the history matching process. If saturation or pressure is changing without wells activity, the model initializes not under static-equilibrium-condition.

Wrong initialization: run without wells shows saturation change.	Correct initialization: saturation doesn't change during run without wells.

Fluids in place calculation

There is a choice of methods to calculate the initial fluid in place values (the initial oil, water and gas saturations) in each grid block. The choice is selected using the ninth argument of the EQUIL keyword, which specifies the number of sub-layers (N) used to obtain the initial saturations:

Center-point equilibration (N=0)
The simulator sets the fluid saturations in each grid block according to the conditions at its center. This method gives a stable initial state since the phase pressure differences in the simulation are also taken between cell centers. It is however the least accurate method, particularly in cases where the fluid contact passes through large grid blocks.

Horizontal fine grid equilibration (N<0)
The top and bottom halves of each grid block are each divided into -N layers of equal thickness, and the saturations are determined locally in each layer. The phase saturations for the block are the average of the saturations in each layer. This option provides a more accurate calculation of the fluids in place, but may yield a solution that is not completely stable in the initial state.

Tilted block fine grid equilibration (N>0)
This option is similar to the previous option, but it takes into account the slope of each grid block. The top and bottom faces of the blocks are treated as planes that are tilted about their central points. The top and bottom halves of the grid block are each divided into N horizontal layers of equal thickness, but the thickness of the layers in the top half is generally different from the thickness of the layers in the bottom half. This is because the distance from the block center line to the highest corner is not the same as the distance from the block center line to the lowest corner, if the upper and lower faces have different tilts. The phase saturations in each block are calculated as a weighted average of the saturations in the layers, weighted according to the area of each layer that is enclosed within the block multiplied by the layer thickness. This option provides the most accurate calculation of fluids in place, but may yield a solution which is not completely stable in the initial state.

With fine grid equilibration there is a redistribution of fluids between grid blocks near the contacts when the simulation begins, which occurs independently of any external driving force (wells etc.). The reason is that a steady state solution on the fine equilibration grid (in which each block is subdivided into several layers) is not necessarily a steady state solution on the coarser simulation grid.

The redistribution of fluids can produce a significant transient when the simulation is started. In ECLIPSE 100 you can overcome this by setting the „quiescence switch' (the switch QUIESC in keyword EQLOPTS). If this switch is on, modifications are applied to the phase pressures to make the initial solution a true steady state. These pressure modifications are applied for the duration of the run.

If EQUIL(9) is not zero, we may therefore have an unstable initial solution, so if the model with no wells and no aquifers is not stable:

Check the 9th argument of the EQUIL keyword
Make sure EQUIL(9) = 0 and run again.

Alternatively in ECLIPSE 100 use QUIECENCE and run again.

Still unstable?
Are you using SWATINIT?

If the initial saturations defined in SWATINIT cannot be honoured then the model may not be stable
This can happen if you are using PPCWMAX with a large initial Sw a long way above the owc

To check for aquifer stability:

Run the model with aquifers but with no wells
Use default values for the initial aquifer pressure
Start with small timesteps and increase the timestep size
Check that there are no convergence problems
Check that field pressures and fluids-in-place are stable Plot SUMMARY quantities such as FPR, FOIP, …
Do not worry about minor changes in FPR, FOIP, ….

SUMMARY section

Never use RPTONLY in the SUMMARY section if you are trying to find the cause of a problem. RPTONLY will only output summary information at report steps and not at timesteps. You may miss some important changes between timesteps.
Add PERFORMANCE to the SUMMARY section, and plot cpu time, elapsed time, timestep size, number of non-linears and number of linears. This can help you find the time during the simulation when the problems are happening, and may help you find the reason.

SCHEDULE data

Multi-Segment Wells

Multi-Segment Wells are more difficult to solve than conventional wells.

Add WSEGITER, with defaults
If you have used TUNING (or any tuning keyword), make sure WSEGITER is after the TUNING keyword.
WVFPEXP – If using VFP tables, can help for THP controlled wells especially with fluctuating gas production.
WELSPECS(10) – turn off cross-flow.
WELSPECS(12) – switching to AVG may reduce instability caused by explicit solution
WELSEGS(8) in record 1 – if using drift flux model, try:

(a) switching to Homogeneous flow
(b) setting OPTIONS(77)

WELSEGS(8) in record 2 – check if friction factor is too high. Lowering this may help convergence.
OPTIONS(30) – governs methods aimed at improving convergence by chopping the increment in all the well‟s solution variables.
OPTIONS(63) – can help if explicit calculation of the hydrostatic head between the well connections and the corresponding grid block centers causes oscillations

Badly/incorrectly defined Multi-Segment Wells definition can be the root cause of extreme convergence problems or even crashes!
Defining MSWs by hand is extremely error-prone. Would only attempt for perfectly horizontal or vertical wells passing through regular box shaped grids. However, best to avoid doing by hand and use Petrel RE instead.

Don‟t do it by hand
If you do, don‟t rush it:

Do it in stages: run the model with standard wells (no segments) to make sure simulation converges, then add the MSW definition.
Draw the segments out on paper on the grid with exact numbers for depth and length before defining the MSWs
Allow 1-2 days to define each well, more if it is complex.
May be easier to do incremental definition than entering absolute numbers when defining the segments by hand. (WELSEGS item 5 of record 2)
To have a segment per cell, need to have n segments, where n is the number of grid cells through which the well will pass. Segment 1 is the zero-tubing head segment.

VFP tables

Use EXTRAPMS in ECLIPSE BLACKOIL. In all simulators, avoid VFP extrapolations.

Other

Does the model run better with no tuning keywords?
Remove TUNING or TSCRIT or CVCRIT keywords and try to run.

Radial runs:
Is the inner radius too small?

Dual porosity runs:
Is sigma too large?