LBL logo Lawrence Berkeley National Laboratory ยป Earth Sciences Division
Questions & Comments | Privacy & Security Notice

TOUGH2 V2.0: Bugs and Fixes

The following list shows known bugs in V2.0 of TOUGH2 and describes how they can be fixed. Holders of a TOUGH2 V2.0 license are strongly encouraged to implement all these fixes into the appropriate source code files and recompile the simulator. All these bug fixes have been implemented in V2.1 of TOUGH2, distributed by the Berkeley Lab Software Center. This list is therefore no longer updated.

go to TOUGH2 home page

June 1, 2007
TOUGH2/ECO2N

Bug. When engaging permeability change due to precipitation or dissolution of salt by specifying IE(11) > 0, the code calculates and reports permeability change (parameter "k-red" in the output file), but the change is not actually implemented in the flow calculations.

June 1, 2007

Fix. The problem is caused by a mismatch of names in subroutine MULTI of TOUGH2 module t2f.f. To correct the problem, change the statement

      if(eosn(1).eq.'EWASG     ') then

in line 2911 to

      if(eosn(1).eq.'EWASG     '.or.eosn(1).eq.'ECO2N     ') then

 

April 27, 2007
module eos9.f

Bug. In eos9.f, two parameters used for internal generation of gravity-capillary equilibrium initial conditions were not initialized when no 'REFCO' specifications were used.

April 27, 2007

Fix. In the DO 33 N=1,NM loop near the top of subroutine EOS, add two lines of coding to initialize parameters ZREF and MRP, as follows.

        DO 33 N=1,NM
        if(mat(n).eq.'REFCO') then
        ...
        ...
        else
           ZREF=0.
           MRP=0
        endif

 

March 9, 2006

CO2TAB fluid property table for module eco2n.f; routine co2tab3.f for generating CO2 property tables

Bug. For low temperatures (< 25 deg-C), the density, viscosity, and specific enthalpy data in the CO2TAB fluid property table supplied with the ECO2N module had lower accuracy at the very highest pressure point (600 bar). Errors reached 3.5 % for density, almost 10 % for viscosity, and 1.5 % for specific enthalpy.

The co2tab3.f routine that was supplied with the ECO2N package for generating fluid property tables has deteriorating accuracy at low temperatures, when pressures approach or exceed 600 bar.

March 9, 2006

Fix. The deterioration in accuracy arose from a flaw in the search interval for the iteration used in co2tab3.f to find the real gas compressibility factor Z. Click here for downloading an improved routine co2tab3.f for generating CO2TAB, in which the search interval was revised to avoid this problem. The revised co2tab3.f follows the iterative bisection process for determining Z with a Newtonian iteration.

Click here to download the updated CO2TAB file generated with the revised co2tab3.f. Please let us know of any problems by e-mail to K_Pruess@lbl.gov

Note that the updated CO2TAB has small differences from the previous version throughout, even for the table points at lower pressures that had satisfactory accuracy. The differences arise because the iteration process for Z is now slightly different. This will have small effects on simulation results, but it may lead to some changes in the behavior of a simulation run, where occasionally convergence now is achieved in a different number of iterations, so that subsequently the automatic time step control may lead to a different choice of delta-t. Once that happens, time truncation errors will of course be different. For testing proper installation of ECO2N by comparing results for sample problems with those given in the user's guide, we recommend using the original CO2TAB.

 

March 8, 2006
module t2f.f

Bug. Occasionally severe convergence problems are encountered when using the MINC method to simulate multiphase flow problems in fractured reservoirs.

March 8, 2006

Fix. The problem occurs because of a conditional switching of relative permeabilities at interfaces in which a nodal distance is zero, as is the case for fracture nodal distances generated by the MINC method. It is recommended to simply comment out (or delete outright) the following two lines in subroutine MULTI, module t2f.f:

        IF(D1.EQ.0..AND.REL20.NE.0.) REL1=REL2
        IF(D2.EQ.0..AND.REL10.NE.0.) REL2=REL1

(These two lines occur a few lines after the following line:)

C-----PERFORM APPROPRIATE UPSTREAM WEIGHTING FOR MOBILITIES.

Note that D=0 triggers other switches as well that should not be changed. These other switches are "static," depending only on input parameters, not on parameters such as relative permeabilities that may change during a simulation.

 

July 8, 2004  Fluid injection at time-dependent rates.

Bug.
When specifying injection of a fluid component at time-dependent rates (LTAB>1 in GENER), and only giving tabular data for times and rates but not for specific enthalpies (ITAB in record GENER.1 left blank), no injection occurs even though non-zero injection rates will be reported in printed output.

July 8, 2004

Fix.
When using time-dependent injection, it is mandatory that tabular data for specific enthalpy of injected fluid be provided as well (array F3, record GENER.1.3). ITAB in record GENER.1 must be specified as a non-blank character in this case, in order that these enthalpy data will be read. Note that when a fluid component is injected at "small" rates (tracer), specific enthalpies may for convenience be specified as zero.

 

January 9, 2002
module t2voc.f

Bug. Subroutine DIFFU did not have a "return" statement, causing bugs on some compilers.

January 9, 2002

Fix. Insert a "return" statement in front of "END" in subroutine DIFFU.

 

September 7, 2001
module ewasg.f

Bug. The calculation of internal energies in BALLA uses the variable NLOC which is undefined. (This will only affect the reporting of global energy balances, not simulation results for grid blocks.)

September 7, 2001

Fix. In subroutine BALLA, module ewasg.f, near the top of the DO 10-loop, insert a statement

        nloc=(n-1)*nk1

in front of

        NLOC2=(N-1)*NSEC*NEQ1

 

April 26, 2001
module ewasg.f

Bugs. In single-phase gas conditions,

(1) vapor pressure is undefined when no non-condensible gas is present, and

(2) gas phase enthalpy calculation has a divide by zero when no vapor is present.

April 26, 2001

Fixes.

(1) Before line 1061 (CALL SUPST(TX,PX,DS,US)) insert the following line

             ps=px

 

 

(2) Replace line 1072 (PAR(NLK2+5)=X3*HAIR+(1.E0-X3)*(US+PS/DS)) with the following two lines

          PAR(NLK2+5)=X3*HAIR
          if(ds.gt.0.) PAR(NLK2+5)=PAR(NLK2+5)+(1.E0-X3)*(US+PS/DS)

 

April 23, 2001
module t2f.f

Bug. Affects EWASG-use only. On November 15, 2000 a revision was introduced whereby permeability change due to solid precipitation is now transmitted through location PAR(...+3) in the solid phase section of the PAR-array, instead of the previously used location PAR(...+6). In order to work properly this revision must also be implemented in subroutines PHASD, PHASF, and GCOR of module t2f.f, which had not been mentioned in the original suggestion of November 15.

April 23, 2001

Fix. In module t2f.f, in lines 3718 (PHASD), 3836 (PHASF), and 4202 (GCOR), change coding from PAR(...+6) to PAR(...+3).

 

January 8, 2001
module t2f.f

Bug. When using fully-coupled multiphase diffusion (default; MOP(24)=0) with the EWASG fluid property module, diffusive fluxes are erroneously printed as 0. This error affects only the printing of fluxes; the simulation uses proper fluxes, and the results calculated for mass fractions etc. are all o.k.

January 8, 2001

Fix. Subroutine OUTDF in file t2f.f has two "PRINT 5071 ..." statements, the first of which is engaged for MOP(24)=1, while the second is active for MOP(24)=0. In the second statement, change the argument of array FDIF from II(KK,NPH,N) to II(KK,2,N).

 

January 3, 2001
module t2cg2.f

Bug. There have been a number of instances where the preconditioned conjugate gradient package T2CG2 gave erratic or unreliable results, occasionally failing altogether. Such failures have occurred when the Jacobian matrix had a large and varying number of zero diagonal elements, and when a grid block was involved in many connections.

January 3, 2001

Fix. A revised version of T2CG2 is available for downloading. Tests of this new version have shown satisfactory performance for problems where the previous version failed, but tests so far have been far from exhaustive. Click here to download a beta-testing version of the revised module (file t2cg22.f, which replaces t2cg2.f of the TOUGH2 V 2.0 package), and a revised INCLUDE file to go with this module (file T2). Please let us know of any problems by e-mail to K_Pruess@lbl.gov

 

November 15, 2000
module ewasg.f

Bug.
Using permeability modification (see section 7.5 in the TOUGH2 User's Guide, LBNL-43134) with EWASG will cause errors in the computation of permeability change from solid precipitation.

November 15, 2000 (amended April 23, 2001)

Fix.
The problem arises from using location PAR(...+6) in the solid phase section of the PAR-array for transmitting permeability change. This should be changed to PAR(...+3), which requires some small modifications in subroutine MULTI, PHASD, PHASF, and GCOR (module t2f.f), and in EOS and OUT (module ewasg.f). Click here to see a text file that summarizes the modifications needed.

 

October 19, 2000
module t2f.f

Bug.
Tabular data for time-dependent heat generation were ignored in Version 2.0; time-dependent mass generation data worked fine.

October 19, 2000

Fix.
The capability to use time-dependent heat generation data can be restored by adding a few lines of code in the HEAT-section of subroutine QU (file t2f.f). Below we show the modified section; the additional lines inserted to enable tabular HEAT generation data are shown bold-faced:

    2 continue
C     ************************************************************
C     * come here when withdrawal or injection of heat           *
C     * at prescribed rate is specified                          *
C     ************************************************************
c
c
C-----IGNORE HEAT SINKS/SOURCES WHEN NOT SOLVING ENERGY EQUATION.
      if(neq.lt.nk1) goto 3900
c
      if(ltaba.gt.1) then
C     ************************************************************
C     * come here when a table of time-dependent generation      *
C     * data is specified (LTABA > 1)                            *
C     ************************************************************
c
         call ttab(n,ltaba,itaba,gn,egn)
         if(igood.ne.0) return
         g(n)=gn
      endif
c

         gpo(n)=g(n)
         r(jmn)=r(jmn)-fac*g(n)
      goto 3900

No other changes are needed. Note that three alternative interpolation schemes may be selected for tabular data through the MOP(12) parameter; see p. 164 in the user's guide (LBNL-43134).

It is suggested that, when making this revision, users should also update the self-documenting date information near the top of QU. The line with the most recent date should be commented out, and a line should be added showing October 18, 2000, as last modification date, as follows.

c 899 FORMAT(/6X,'QU       1.1      23 January   1998',6X,
  899 FORMAT(/6X,'QU       1.1      18 October   2000',6X,

 

July 6, 2000  T2VOC

Bug.
For T2VOC sample problem *rgdif*, vapor diffusivity (printed as "DIFFW") is slightly different in Version 2.0 as compared to the previous version of T2VOC.

July 6, 2000

Fix.
The discrepancy was traced to the fact that in the earlier version of T2VOC there is a default value of 1.8 for the temperature exponent TEXP of diffusivity (record PARAM.1), while in the T2VOC that is part of TOUGH2 V2.0 the actual input value is used, even if this value happens to be zero. No coding change was made; users just need to be aware of this change in default settings.

 

May 25, 2000  T2VOC

Bug.
The parameter "FOC" (fraction organic carbon) read in record ROCKS.1.1 was not properly initialized in TOUGH2, V 2.0.

May 25, 2000

Fix.
Insert a statement "foc(nm)=xkd3(nm)" after line 210 in module t2f.f, as shown in bold face, below:

209       IF(NAD.GE.1) READ 7,COM(NM),EXPAN(NM),CDRY(NM),TORT(NM),GK(NM)
210      x,XKD3(NM),XKD4(NM)
211 c.....5-25-00: assign fraction organic carbon for T2VOC
212       foc(nm)=xkd3(nm)
213     7 FORMAT(8E10.4)