The calculated secondary variables are used to assemble the linear equations (Eq. 5-9) that are solved at each step of the Newton-Raphson iteration. The structure of the coupled linear equations at iteration step (p+1) is shown in Figure 13. When a time step is converged, the primary variables X are updated, X X + DX. These equations are arranged and numbered sequentially, with the first NK equations per grid block representing component mass balances, while the last equation represents the energy balance. The row indices of the Jacobian matrix correspond to the component balance equations, while the column indices correspond to the sequence of primary variables in array X.
In many applications, heat effects and temperature changes may be small, so that it would be sufficient to consider just the mass balance equations. To force temperature changes to zero, we can assign the solid matrix an overwhelmingly large heat capacity, so that finite rates of heat exchange will cause negligible temperature change. This approach is perfectly acceptable and useable with TOUGH4, but it has the disadvantage that the full number of balance equations must be solved, even though the energy balance reduces to the trivial statement ∂T/∂t = 0. TOUGH4 offers an option to run problems without temperature changes, at considerable savings of computing time. TOUGH4 uses a parameter NEQ, distinct from NK, to count the number of balance equations to be solved per grid block. Normally we would have NEQ = NK + 1. However, for certain EOS modules that have temperature as primary variable # NK1, users can choose to assign NEQ = NK in the input file; then no energy equations are set up or solved and the number of coupled equations is reduced by one per grid block. In this case, only the first NK primary variables per grid block will contribute matrix columns, while variable # NK1, which must be temperature, remains passive and is not engaged or altered in the linear equation handling. However, all thermophysical parameters will be calculated at the temperature values specified in variable # NK1. Difference to the previous version of TOUGH, TOUGH4 allows flexible change of NK at runtime based on the number of gas components, tracers/salts, and whether including second water/brine in current simulation. This may achieve additional efficiency by solving fewer equations.
Note that the accumulation terms of the balance equations depend only on primary variables for one grid block, so that they will generate non-zero derivative terms only in an NEQ NEQ submatrix that is located on the diagonal of the Jacobian matrix J. The flow terms, being dependent on primary variables of two grid blocks, will generate two non-zero NEQ NEQ submatrices of derivatives, which are located in the off-diagonal matrix locations corresponding to the two grid blocks.
In TOUGH4 all Jacobian matrix elements as well as the entries in the vector rhs of residuals are calculated in subroutine Jacobian_SetUp. The calculation first assigns all matrix elements arising from the accumulation terms, of which there are NEQ NEQ. Calculation of the derivatives demands that each accumulation term is calculated NEQ + 1 times; once for the state point (X + DX), and NEQ times for each of the NEQ primary variables incremented (X + DX + DELX). Additional contributions to diagonal terms in the Jacobian J may arise from sink and source terms, if present; these are assigned in subroutine WellFlux. The sink and source terms also include contribution from other processes, such as biodegradation reaction (by subroutine Bio_react) and heat exchange from analytical solution (by subroutine AnalyticalHeatExchange). Subsequently all flux terms are evaluated. These depend in general on the 2·NEQ primary variables of the two connected grid blocks, so that a total of 2·NEQ + 1 flux terms need to be evaluated for calculation of the state point as well as of all derivative terms.
After all matrix elements and members of the right-hand side vector of residuals have been assembled, the linear equations (Eq. (5-9)) are solved by user-selected linear solver. The resulting increments in the primary variables are added to the array DX, and the process of linear equation setup and solution is repeated for the primary variables X + DX. This process continues until the residuals are reduced below a preset convergence tolerance (Eq. (5-10)). If convergence is not achieved within a specified maximum number of iterations (usually 8, see NOITE in ), the time step is repeated with reduced time increment
