Several STOMP Operational Modes couple energy and mass transport. Energy transfer is similar to mass component transfer, occurring through advection and diffusion of mass as well as thermal diffusion through the fluid and solid phases.
Table 1. STOMP operational modes that consider heat transfer in liquid (l), gas (g), NAPL (n), and solid (s) phases.
The energy conservation equation equates the time rate of change of energy within a control volume with the flux of energy crossing the control volume surface and energy sources associated with mass or heat sources.
Where the energy accumulation term is
and heat transfer can occur in the aqueous (l), NAPL (n), gas (g), hydrate (h), ice (i), solid (s), and precipitated (p) phases under equilibrium conditions, depending on the operational mode.
Thermal flux is a combination of advective and diffusive (or conductive) components.
Surface integrals are approximated by discretizing the control volume surfaces into node surfaces and summing the contributions to thermal flux over the node surfaces. The energy accumulatoin term (i.e., lefthandside terms) are summed and integrated over the node volume according to
Defining flux directions parallel to the surface normal allows the surface integrals to be converted to summations over all node surfaces.
This transformation strictly requires an orthogonal grid system for the flux directions to be aligned with the surface normals. Nonorthogonal systems will yield energy balance errors.
The mass conservation equations are discretized in time using a fully implicit scheme, where the time levels are indicated with superscripts. The primary unknowns for the mass conservation equations are intrinsic properties at node volume centroids (node grid point) for time level t+δt.
The residual equation for energy is then the difference between lefthandside and righthandside.
In the STOMP simulator, thermal energy is partitioned, according to thermal equilibrium conditions, among the fluid and solid phases. The thermal capacitance of unconnected pore space, represented by the difference between the total and diffusive porosity, is computed as it is filled with liquid water.
Heat transfer by hydraulic dispersion of flowing fluid phases is neglected.
Enhanced vapor transport is incorporated through enhancement factors for component diffusion through the gas phase.
Energy associated with component mass sources are included as internal heat generation sources.
Reference states for enthalpy and internal energy are component dependent.
Latent heat transport is considered through vapor transport through the gas phase and equilibrium thermodynamics.
STOMPHYDTKE ignores energy flux associated with component diffusive fluxes.
Symbols(In order of appearance) 

time, s  
volume of element n, m^{3} 

energy accumulation term for component j , J/m^{3}  
surface of element n, m^{2}  
heat flux tensor, W/m^{2}  
unit surface normal vector  
specific enthalpy of phase γ, J/kg  
specific mass source of phase γ, kg/m^{3} s  
specific heat source, W/m^{3}  
diffusive porosity 

density of phase γ, kg/m^{3} 

saturation of phase γ  
specific internal energy of phase γ, J/kg  
total porosity  
effective thermal conductivity, W/Km  
temperature, K  
Darcy velocity vector of phase γ, m/s  
diffusivedispersive flux of component j for the phase γ, kg/m^{2} s  
relative permeability of phase γ  
intrinsic permeability, m^{2}  
kinematic viscosity of phase γ, Pa s  
pressure of phase γ, Pa  
acceleration of gravity, m/s^{2}  
unit gravitational direction vector  
molecular weight of component j, kg/kg mol  
molecular weight of phase γ, kg/kg mol  
phase tortuosity for phase γ  
diffusion coefficient of component j for phase, m^{2}/s  
mole fraction of component j in phase γ  
area of surface, m^{2} 
Subscripts  

phase index  
aqueous liquid phase  
nonaqueous liquid phase  
gaseous phase  
hydrate phase  
ice phase  
precipitated salt phase  
solid phase  
node surface index  
negative surface (west, south, bottom)  
positive surface (west, south, bottom) 
Superscripts  

component index  
upwind weighting scheme  
harmonic weighting scheme 