I will share my thoughts and tips on reservoir simulation (TOUGH+RealGasBrine, TOUGH2, TOUGHREACT, TOUGH+HYDRATE) here from time to time.
Welcome to discuss with me, or contact me if you would like to post your thoughts here.

Nonlinear Convergence Problems

Under challenging conditions (for example, highly nonlinear behavior associated with phase appearance and disappearance), residual reduction may be difficult and convergence may not be attained within a reasonable number of iterations. This often triggers timestep cutbacks to impractically small levels. Convergence problems can be categorized as mathematical or physical. Given the robustness of TOUGH+, most issues stem from inappropriate physical settings. A clear understanding of the physical process before starting simulations cannot be overemphasized.

  1. My test results indicate that the stabilized BiConjugate Gradient (BCGS) solver combined with the block Jacobi (BJACOBI) preconditioner consistently provides the most robust and efficient performance across elliptic, parabolic, and hyperbolic governing equations, and is therefore recommended as the default configuration.
  2. Phase changes are the most common reason for nonconvergence. If the system is very sensitive to phase change but the variation in primary variables is not large enough to support a complete transition, local phase conditions can oscillate between two thermodynamic states, resulting in very small timesteps. TOUGH+ introduces overshoot variables to trigger phase and state changes more robustly. Proper settings allow one-way phase changes instead of repeated oscillations.
  3. Phase changes, if improperly described, can induce numerical discontinuities that prevent Newton-Raphson convergence. For example, relative permeability or capillary pressure curves with steep gradients and no endpoint smoothing can create derivatives that exceed compiler limits, causing the Jacobian to become singular when a phase appears or disappears. Ensure derivative continuity by using appropriate smoothing that remains faithful to the data.
  4. Incorrect derivatives inevitably lead to successive timestep cutbacks to minuscule levels before the simulation is suspended automatically. In practice, repeated timestep reductions are a clear indicator of wrong derivatives.
  5. Simulations can appear to struggle with convergence as the system approaches steady state. In this case, Jacobian derivatives tend toward zero, and the solution is effectively steady with minimal changes in primary variables. This is not a convergence failure, but an indication of steady state, often accompanied by rapidly increasing timesteps.
  6. It is easy to identify steady-state behavior for the whole system, but it can be challenging to detect it locally. For example, when gas plumes migrate upward, grid blocks at the lateral edges can constrain the timestep because inflow and outflow fluxes nearly balance. These blocks struggle with phase transitions, so a relaxed phase-change criterion may be appropriate if local transitions are not critical to global results.
  7. Set the global convergence criterion appropriately, especially for large-scale simulations with elements of vastly different volumes. If it is too small, the model requires very small timesteps and long runtimes; if too large, accuracy may suffer. The default relative criterion is 1e-5. For very large simulations, 1e-4 may be acceptable, but larger values are discouraged. A looser criterion is preferable when the model has (i) millions of degrees of freedom, (ii) long simulation times (for example, CO2 storage over 50 years), and (iii) smooth or quasi-steady solutions without sharp gradients.
  8. [TOUGHREACT, personal communication with Dr. Eric Sonnenthal]
    • Capturing natural system behavior improves convergence. Deciding which minerals and species to remove depends on model purpose (for permeability changes only, or production chemistry as well). A species with low concentration can still be important if it has low solubility and quickly reincorporates into secondary minerals.
    • Removing a few less-important minerals usually has little effect; convergence issues typically depend on the more important minerals.
    • Many convergence issues come from placing far-from-equilibrium fluids into hot rock instantaneously. For high-temperature problems, run a 1-block pseudo-steady-state simulation starting at low temperature and step up in ~50C intervals to reservoir temperature. Use the resulting fluid composition as input to the large model. This initialization is critical. Progress carefully from TH to THM to THMC to obtain native-state conditions before injection.
    • Grid refinement does not strongly affect chemical equilibration, but it does affect time-stepping through fluid velocity. Time steps should be constrained by the Courant limit and reaction rates, which are usually much smaller than typical fluid-flow models.
    • Reactive surface areas are extremely important and should be calculated based on the conceptual model, then calibrated to native-state conditions.
    • If the Courant limit is much greater than 1, results can be erroneous. Reactive chemistry is not like fluid flow or stress calculations; small concentration errors (for example, 1e-4 or 1e-5) accumulate and directly affect porosity and permeability. If the timestep exceeds the Courant limit, concentrations become unreasonable, convergence degrades, and the run can crash.

Mesh Generation

Mesh generation is critical in numerical simulation because it directly affects accuracy, stability, and efficiency. Mesh quality determines how well geometry, boundary conditions, and physics are captured. A well-constructed mesh resolves key regions with high gradients or complex behavior, avoids large aspect ratios, and prevents abrupt transitions between coarse and fine elements.

  1. MESHMAKER V1.5 (Moridis 2016) can generate two-dimensional radially symmetric (r, z) meshes and one-, two-, and three-dimensional rectilinear (Cartesian) grids. It produces grids directly usable by TOUGH+ and TOUGH2 family codes.
  2. TOUGHIO (Luu 2020) is an open-source library that provides tools for TOUGH pre- and post-processing using Python. It can import meshes from external software (Abaqus, FLAC3D, Gmsh) and convert them to TOUGH MESH files. This enables flexible mesh preparation. In this study, some meshes are created in Gmsh (Geuzaine and Remacle 2009) and converted to TOUGH+ MESH files. Note that a good mesh must satisfy orthogonality conditions, so quality checks are essential, especially for unstructured meshes.

    The recommended workflow is: (i) design the mesh in Gmsh and export *.msh and *.vtk files; (ii) check mesh quality using the pyvista library in TOUGHIO; (iii) convert the *.msh file to a TOUGH MESH file using TOUGHIO read and write functions; and (iv) define system boundaries in the MESH file.

Compiler

  1. [TOUGH+] When simulating Discrete Fracture Models with very small element volumes (for example, 1.0e-8 in fracture cells), compiling with the Intel Compiler option '-O3' can cause convergence issues. The more aggressive optimizations can affect floating-point precision. Consider compiling without optimizations for such cases. I have not confirmed whether the same issue occurs with GNU.
  2. [TOUGH+] Release builds may not initialize variables, leading to unpredictable behavior. Add options to initialize variables to zero for safety.
  3. [TOUGH+] Using an initial value of '0' in the INCON file is risky and often causes convergence issues, especially in large models. Use a very small value such as 1.0e-18 instead.
  4. When using parallel code on Linux, set GFORTRAN_UNBUFFERED_ALL=1 if you need unbuffered output for immediate printing.
  5. [TOUGHREACT] If a segmentation fault occurs immediately after launch, increase the OMP_STACKSIZE environment variable or reduce the number of elements and connections in flowpar_v4.inc to lower memory requirements.