User login

You are here

Journal Club Theme of September 2014: Numerical modeling of thermo-hydro-mechanical coupling processes in porous media

WaiChing Sun's picture

Thermo-hydro-mechanics (THM) is a branch of mechanics aimed to predict how deformable porous media behave, while heat transfer and fluid transport simultaneously occur in the pores filled by liquid and/or gas. Understanding these multi-physical responses is important for a wide spectrum of modern engineering applications, such as tissue scaffolding, geothermal heating, mineral exploration and mining, hydraulic fracture, energy piles, tunneling with frozen soil and nuclear waste storage and management. The major difficulty encountered for modeling this coupling physical processes is that the thermal and pore-fluid diffusions, and the deformation of solid skeleton may take place in profoundly different spatial (micron vs. kilometer) and temporal scales (second vs decades). Bundling all the physical processes in simulations while maintaining accuracy and robustness is therefore a difficult task.In the last three decades, a considerable progress has been made for deriving mathematical theories and implementing computer models to replicate the fully coupled thermo-hydro-mechanical processes. This journal club article will focus on the derivation and implementation of the thermo-hydro-mechanical model for numerical simulations. Any feedback or critique are greatly appreciated.

1. Monolithic vs. operator-splitting schemes

The governing equation of the thermo-hydro-mechanical problem can be obtained from the balance principles. The displacement, pore pressure, Darcy’s velocity and temperature are usually the primary solution field stored at the nodal points. If the whole set of governing equations are solved simultaneously, then the solver is called monolithic. The monolithic solver is particularly important for strongly coupling problems in which passing parameters among fluid and solid solvers may not be sufficient (Preisig & Prevost, 2011).For instance, a monolithic small strain finite element code, FRACON, has been developed by (Nguyen & Selvadurai,1995). In this code, the balance of linear momentum and mass are fully coupled, while thermal transport may affect the solid deformation and pore-fluid diffusion but not vice versa. Li, Liu and Lewis, introduce a co-rotational FEM formulation and incorporate plasticity into THM model to model the non-isothermal elastoplastic responses of porous media at large strain in (Li,2005). Recent work by Preisig and Prevost employed a fully coupled THM simulator to compare the numerical solutions against the field data in a case study for carbon dioxide injection at In Salah, Algeria (Preisig & Previst, 2011). Kolditz et al., 2012 introduces an open-source project OpenGeoSys, which takes advantage of an object-orient framework and provides software engineering tools such as platform-independent compiling and automated benchmarking for developers.

In addition to the monolithic finite element scheme, attempts have been made to sequentially coupled multiphase flow and geomechanical simulators by establishing proper feedback and information exchange mechanisms. One such example is TOUGH-FLAC, which links flow simulator TOUGH2 with a small strain finite difference code FLAC (Rutqvist,2011). Klar et al. 2013 employ an explicit time-marching scheme to simulate the fully coupled thermal, and pore-fluid diffusions and the path dependent responses of the solid skeleton. This sequential coupling approach is an attractive alternative to the monolithic approach, as it is easier to implement and maintain flow and solid simulators separately. In addition, the operator-splitting approach, if used properly, can offer tremendous saving on memory and CPU times. The operator splitting approach also enables one to assign different time steps to different physical processes. For instance, in reactive-diffusion simulations, it is common to have the reaction component updating with a several order finer time step than the diffusion counterpart. However, proper communication must be established to ensure the correctness and numerical stability of numerical solutions (Jha & Juanes, 2007, Preisig, 2011, Klar et al. 2013, Sun et al., 2013a, Sun et al. 2013b). As noted in (Jha & Juanes,2007), the sequential coupling scheme used to link the fluid and solid simulators may have profound impact on the efficiency,  stability and accuracy of the numerical solutions (Schrefler et al.,1995, Schrefler et al.,1997).

If the fluid and solid simulators use different grids or meshes (eg. finite volume for fluid and finite element for solid), then a proper data projection scheme is required to transfer information from Gauss points and nodes of the solid mesh to the fluid mesh and vice versa (Goumiri & Previst, 2011). For large scale parallel simulations, the sequential couplings must be carefully designed to avoid causing bottleneck due to the difference in solver speed.  This can be a significant problem if either the solid or fluid solver runs only in serial. 

2.  Mixed finite elements (e.g. Taylor-Hood, Raviart-Thomas)  vs.  stabilized procedures

As noted in  Zienkiewicz et al., 1999, Wan, 2002, Mira et al., 2003, Truty & Zimmermann, 2006, Simoni et al., 2008, White & Borja, 2008, Preisig & Prevost , 2011, Sun et al, 2013a, Sun et al. 2013b Borja, 2013, numerical stability is often a major challenge for poromechanics models. Due to the lack of inf-sup condition (Babuvska, 1973, Brezzi, 1985, Bathe, 2001, Bochev et al. 2006),  pore pressure and temperature fields may exhibit spurious oscillation patterns and/or checkerboard modes if the displacement, pore pressure and temperature are spanned by the same set of basis function. While these spurious oscillations are less severe at the drained/isothermal limit, they may intensify when small time step is used or when materials are near undrained/adiabatic limit. From a mathematical viewpoint, these non-physical results are due to the kernel (null space) of the discrete gradient operator being non-trivial. This non-trivial kernel makes it possible to have certain spatially oscillating pore pressure and temperature fields that have no impact on the solid deformation in a numerical simulation.

To cure the numerical instability due to the lack of inf-sup condition, previous researches have established a number of techniques that employ different basis functions to interpolate displacement and pore pressure and obtain stable solutions.  For instance, Zienkiewicz and coworkers Zienkiewicz et al. 1999, and Borja  & Alarcon 1995 used Taylor-Hood finite element (quadratic basis functions for displacement and linear basis functions for pore pressure) to satisfy inf-sup condition and maintain numerical stability for isothermal hydromechanics problems.

On the other hand, Jha & Juanes 2007 have shown that linear displacement combined with pore fluid velocity in the lowest-order Raviart-Thomas space, and piecewise constant pore pressure may also lead to stable solutions for isothermal poromechanics problems. Nevertheless, inf-sup stable mixed finite element models require multiple meshes for displacement, pore pressure and/or fluid velocity. As a result, they require additional programming effort to pre- and post-processing data and maintain the more complex data structure.

To avoid the complications of using multiple meshes for each solution field, an alternative is to use one single equal-order finite element mesh with stabilization procedures. Many stabilization procedures have been proven to be able to eliminate the spurious oscillation modes without introducing extra diffusion for small strain isothermal poromechanics problems. For instance, White & Borja 2008 employed a polynomial projection scheme originated from Dohrmann & Bochev, 2004 to simulate slip weakening of a fault segment. A numerical example in this paper shows that spurious oscillation may persist in the inf-sup stable finite elements (e.g. quadratic-displacement/linear-pore-pressure pair) if the diffusivity is very low for a given time step size. Perhaps the major drawback for the stabilized finite element approach is that it is often necessary to find the stabilization parameter that just give the right amount of stabilization effect without over-diffusing the solution. Recently, my collaborators and I have attempted to address this issue for isothermal poromechanics problem in the large deformation regime, where the stabilization term is adaptively adjusted to avoid over-diffusion (Sun 2013a and Sun et al. 2013b). The stabilization term is derived by solving a small strain one-dimensional deformation-diffusion problem. By determining the critical value of the diffusivity that leads to complex valued growth/decay rate, the optimal value for the stabilization parameter can be estimated. For the three-field thermo-hydro-mechanical problem, the situation is more complicated. Since the thermal convection-diffusion and pore fluid diffusion may reach steady-state in profoundly different rates, it is unclear how to stabilize both the pore pressure and temperature fields without over-diffusing one or both of them.

3.   Length scale and bifurcations 

Unlike single-phase materials, porous media are multiphase in which solid, liquid(s) and gas(s) can all occupy fractions of volume of the representative elementary volume. As a result, the deformation of the solid skeleton are influenced by the heat transfer and pore fluid diffusions. One direct consequence is that the diffusion will introduce rate dependence to the mechanical responses. This rate dependence is sometime considered by some as a localization limiter which regularizes the governing equations when deformation band formed in numerical simulations. Zhang, et al (1999) analyzed the characteristic equation of a one-dimensional wave propagation problem and found that a length scale may be derived for compressive wave even when softening occurs for solid skeleton, but this is not the case for shear wave. Abellan & de Borst (2005) conducted a similar dispersion analysis and found that the physical length scale disappear in short wavelength limit. This result is supported by numerical simulations in which the strain of a softening bar composed of saturated porous media is found to be mesh dependence. In other words, using diffusion alone as a mean to circumvent mesh dependence is not productive, according to both Zhang et al (1999) and Abellan & de Borst (2005). Nevertheless, a mesh independent post-bifurcation response may still be obtained in numerical simulations, if a length scale can be introduced through other means (e.g. grain rotation, gradient dependence, nonlocal plasticity or damage model, rate dependent solid constituent…etc). 

A similar observation is made by Bazant (2010) for models that couple multiple spatial scales. The author observed that incorporating sub-scale simulations to macroscopic simulations alone does not introduce length scale for the softening materials. This treatment merely move the strain localization phenomenon one scale down. Recent multiscale discrete element/finite element models proposed by Guo & Zhao (2014) and by Nguyen et al. (2014) both confirm that the macroscopic finite element responses are mesh dependence in the post-bifurcation regime even though the macroscopic responses are inferred from grain-scale simulations. 


1.Abellan, M. A., & De Borst, R. (2006). Wave propagation and localisation in a softening two-phase medium. Computer Methods in Applied Mechanics and Engineering, 195(37), 5011-5019.

2.Armero F (1999) Formulation and finite element implementation of a multiplicative model of cou- pled poro-plasticity at finite strains under fully saturated conditions. Computer methods in applied mechanics and engineering 171(3):205–241 

3.Athy LF (1930) Density, porosity, and compaction of sedimentary rocks. AAPG Bulletin 14(1):1–24 

4.Auricchio F, Beir ̃ao da Veiga L, Lovadina C, Reali A (2005) A stability study of some mixed finite elements for large deformation elasticity problems. Computer methods in applied mechanics and 
engineering 194(9):1075–1092 

5.Auricchio F, da Veiga LB, Lovadina C, Reali A, Taylor RL, Wriggers P (2013) Approximation of 
incompressible large deformation elastic problems: some unresolved issues. Computational Mechanics 
pp 1–15 

6.Babuˇska I (1973) The finite element method with lagrangian multipliers. Numerische Mathematik 

7.Bathe KJ (2001) The inf-sup condition and its evaluation for mixed finite element methods. Com- 
puters and Structures 79(2):243 – 252, DOI 10.1016/S0045-7949(00)00123-1 

8.Bazant, Z. P. (2010). Can multiscale-multiphysics methods predict softening damage and structural failure?. International Journal for Multiscale Computational Engineering, 8(1).

9.Bear J (1972) Dynamics of fluids in porous media. Elsevier Publishing Company, New York, NY 

10.Biot M (1941) General theory of three dimensional consolidation. Journal of Applied Physics 
12(2):155 –164 

11.Bochev PB, Dohrmann C, Gunzburger M (2006) Stabilization of low-order mixed finite elements for 
the stokes equations. SIAM J Numer Anal 44(1):82–101 

12.Borja RI (2013) Plasticity Modeling and Computation. Springer, Berlin 

13.Borja RI, Alarco ́n E (1995) A mathematical framework for finite strain elastoplastic consolidation 
part 1: Balance laws, variational formulation, and linearization. Computer Methods in Applied Me- 
chanics and Engineering 122(1):145–171 

14.Borja RI, Tamagnini C, Alarc ́on E (1998) Elastoplastic consolidation at finite strain part 2: Finite 
element implementation and numerical examples. Computer Methods in Applied Mechanics and 
Engineering 159(1):103–122 

15.Brezzi F, Fortin M (1991) Mixed and hybrid finite element methods. Springer-Verlag New York, Inc. 

16.Brezzi F, Douglas J, Marini LD (1985) Two families of mixed finite elements for second order elliptic 
problems. Numerische Mathematik 47:217–235 

15.Broccardo M, Micheloni M, Krysl P (2009) Assumed-deformation gradient finite elements with nodal integration for nearly incompressible large deformation analysis. International journal for numerical methods in engineering 78(9):1113–1134 

16.Callari C, Armero F (2004) Analysis and numerical simulation of strong discontinuities in finite strain poroplasticity. Computer Methods in Applied Mechanics and Engineering 193(27):2941–2986 

17.Castellazzi G, Krysl P (2012) Patch-averaged assumed strain finite elements for stress analysis. International Journal for Numerical Methods in Engineering 90(13):1618–1635 

18.Coussy O (2004) Poromehcanics 

19.Coussy O (2007) Revisiting the constitutive equations of unsaturated porous solids using a lagrangian 
saturation concept. International Journal for Numerical and Analytical Methods in Geomechanics 

20.Cowin SC, Doty SB (2009) Tissue mechanics. Springer 

21.Dohrmann CR, Bochev PB (2004) A stabilized finite element method for the stokes problem based on 
polynomial pressure projections. International Journal for Numerical Methods in Fluids 46(2):183– 

22.Girault V, Raviart PA (1986) Finite element methods for Navier-Stokes equations, Theory and 
algorithms, volume 5 of Springer Series in Computational Mathematics. Springer, Berlin 

23.Goumiri IR, Prevost JH (2011) Cell to node projections: An assessment of error. International Journal for Numerical and Analytical Methods in Geomechanics 35(7):837–845, DOI 10.1002/nag.927, URL 

24.Guo, N., & Zhao, J. (2014). A coupled FEM/DEM approach for hierarchical multiscale modelling of granular media. International Journal for Numerical Methods in Engineering, 99(11), 789-818.

25.Heroux MA, Bartlett RA, Howle VE, Hoekstra RJ, Hu JJ, Kolda TG, Lehoucq RB, Long KR, 
Pawlowski RP, Phipps ET, et al (2005) An overview of the trilinos project. ACM Transactions on 
Mathematical Software (TOMS) 31(3):397–423 

26.Hiroshi H, Minoru T (1986) Equivalent inclusion method for steady state heat conduction in com- 
posites. International Journal of Engineering Science 24(7):1159–1172 

27.Holzapfel GA (2000) Nonlinear solid mechanics: a continuum approach for engineering 

28.Howell JS, Walkington NJ (2011) Inf–sup conditions for twofold saddle point problems. Numerische 
Mathematik 118(4):663–693 

29.Jha B, Juanes R (2007) A locally conservative finite element framework for the simulation of coupled 
flow and reservoir geomechanics. Acta Geotechnica 2:139–153, DOI 10.1007/s11440-007-0033-0, URL 

30.Jing L, Tsang CF, Stephansson O (1995) Decovalexan international co-operative research project on 
mathematical models of coupled thm processes for safety analysis of radioactive waste repositories 

31.Karrech A, Poulet T, Regenauer-Lieb K (2012) Poromechanics of saturated media based on the 
logarithmic finite strain. Mechanics of Materials 51(0):118 – 136 

32.Kolditz O, Bauer S, Bilke L, Bo ̈ttcher N, Delfs J, Fischer T, Go ̈rke U, Kalbacher T, Kosakowski 
G, McDermott C, et al (2012) Opengeosys: an open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (thm/c) processes in porous media. Environmental Earth Sci- ences 67(2):589–599 

33.Lewis R, Majorana C, Schrefler B (1986) A coupled finite element model for the consolidation of nonisothermal elastoplastic porous media. Transport in Porous Media 1(2):155–178 

34.Li X, Liu Z, Lewis RW (2005) Mixed finite element method for coupled thermo-hydro-mechanical process in poro-elasto-plastic media at large strains. International Journal for Numerical Methods in Engineering 64(5):667–708 

35.Liu R, Wheeler M, Dawson C, Dean R (2009) Modeling of convection-dominated thermoporome- chanics problems using incomplete interior penalty galerkin method. Computer Methods in Applied Mechanics and Engineering 198(9):912–919 

36.McTigue D (1986) Thermoelastic response of fluid-saturated porous rock. Journal of Geophysical Research 91(B9):9533–9542 

37.Mira P, Pastor M, Li T, Liu X (2003) A new stabilized enhanced strain element with equal order of in- terpolation for soil consolidation problems. Computer methods in applied mechanics and engineering 192(37):4257–4277 

37.Moran B, Ortiz M, Shih C (1990) Formulation of implicit finite element methods for multiplicative finite deformation plasticity. International Journal for Numerical Methods in Engineering 29(3):483– 514 

38.Nguyen T, Selvadurai A (1995) Coupled thermal-mechanical-hydrological behaviour of sparsely frac- tured rock: Implications for nuclear fuel waste disposal. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 32(5):465 – 479 

39.Nguyen, T. K., Combe, G., Caillerie, D., & Desrues, J. (2014). FEM× DEM modelling of cohesive granular materials: Numerical homogenisation and multi-scale simulations. Acta Geophysica, 62(5), 1109-1126.

40.Notz PK, Pawlowski RP, Sutherland JC (2012) Graph-based software design for managing complexity and enabling concurrency in multiphysics pde software. ACM Transactions on Mathematical Software (TOMS) 39(1):1 

41.Nur A, Byerlee J (1971) An exact effective stress law for elastic deformation of rock with fluids. Journal of Geophysical Research 76(26):6414–6419 

42.Pantuso D, Bathe KJ (1997) On the stability of mixed finite elements in large strain analysis of incompressible solids. Finite elements in analysis and design 28(2):83–104 

43.Pawlowski RP, Phipps ET, Salinger AG (2012) Automating embedded analysis capabilities and man- aging software complexity in multiphysics simulation, part i: Template-based generic programming. Scientific Programming 20(2):197–219 

44.Pawlowski RP, Phipps ET, Salinger AG (2012) Automating embedded analysis capabilities and man- aging software complexity in multiphysics simulation, part i: Template-based generic programming. Scientific Programming 20(2):197–219 

45.Pawlowski RP, Phipps ET, Salinger AG, Owen SJ, Siefert CM, Staten ML (2012) Automating em- bedded analysis capabilities and managing software complexity in multiphysics simulation, part ii: Application to partial differential equations. Scientific Programming 20(3):327–345 

46.Postelnicu A (2004) Influence of a magnetic field on heat and mass transfer by natural convection from vertical surfaces in porous media considering soret and dufour effects. International Journal of Heat and Mass Transfer 47(6):1467–1472 

47.Preisig M, Pr ́evost JH (2011) Coupled multi-phase thermo-poromechanical effects. case study: Co2 injection at in salah, algeria. International Journal of Greenhouse Gas Control 5(4):1055 – 1064, DOI 10.1016/j.ijggc.2010.12.006 

48.Prevost JH (1982) Nonlinear transient phenomena in saturated porous media. Computer Methods in Applied Mechanics and Engineering 30(1):3 – 18 

49.Radovitzky R, Ortiz M (1999) Error estimation and adaptive meshing in strongly nonlinear dynamic problems. Computer Methods in Applied Mechanics and Engineering 172(1):203–240 

50.Rajagopal KR, Tao L (1995) Mechanics of mixtures, vol 754. World scientific Singapore 

51.Rice J, Cleary M (1976) Some basic stress diffusion solutions for fluid-saturated elastic porous media 
with compressible constituents. Rev Geophys 14(2):227–241, DOI 10.1029/RG014i002p00227 

52.Rutqvist J (2011) Status of the tough-flac simulator and recent applications related to coupled fluid 
flow and crustal deformations. Computers and Geosciences 37(6):739 – 750 

53.Salinger AG, Pawlowski RP, Phipps ET, Bartlett RA, Hansen GA, Kalashnikova I, Ostien JT, Sun W, Chen Q, Mota A, Muller RP, Nielsen E, Gao X (2013) Albany: A component-based partial 
differential equation code built on trilinos. ACM Transactions on Mathematical Software 

54.Schrefler B (1995) Fe in environmental engineering: coupled thermo-hydro-mechanical processes in porous media including pollutant transport. Archives of Computational Methods in Engineering 

55.Schrefler B, Simoni L, Turska E (1997) Standard staggered and staggered newton schemes in thermo- 
hydro-mechanical problems. Computer methods in applied mechanics and engineering 144(1):93–109 

56.Selvadurai A, Nguyen T (1997) Scoping analyses of the coupled thermal-hydrological-mechanical behaviour of the rock mass around a nuclear fuel waste repository. Engineering Geology 47(4):379– 

57.Selvadurai A, Suvorov A (2012) Boundary heating of poro-elastic and poro-elasto-plastic spheres. 
Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science 468(2145):2779– 

58.Simo J, Miehe C (1992) Associative coupled thermoplasticity at finite strains: Formulation, numerical 
analysis and implementation. Computer Methods in Applied Mechanics and Engineering 98(1):41– 104 

58.Simo J, Taylor R, Pister K (1985) Variational and projection methods for the volume constraint in finite deformation elasto-plasticity. Computer Methods in Applied Mechanics and Engineering 51(1):177–208 

59.Simoni L, Secchi S, Schrefler B (2008) Numerical difficulties and computational procedures for thermo-hydro-mechanical coupled problems of saturated porous media. Computational Mechanics 43(1):179–189 

60.Skempton A (1984) Effective stress in soils, concrete, and rocks. Selected papers on soil mechanics pp 4–16 Souza Neto EA, Peri ́c D, Owen DRJ (2008) Computational Methods for Plasticity. John Wiley & Sons, Ltd 

62.Sun W, Andrade J, Rudnicki J (2011) Multiscale method for characterization of porous microstruc- tures and their impact on macroscopic effective permeability. International Journal for Numerical Methods in Engineering 88(12):1260–1279 

63.Sun W, Andrade JE, Rudnicki JW, Eichhubl P (2011) Connecting microstructural attributes and permeability from 3d tomographic images of in situ shear-enhanced compaction bands using multi- scale computations. Geophysical Research Letters 38(10):L10,302 

64.Sun W, Chen Q, Ostien J (2013a) Modeling the hydro-mechanical responses of strip and circular punch loadings on water-saturated collapsible geomaterials. Acta Geotechnica pp 1–32, DOI 10.1007/s11440- 013-0276-x 

65.Sun W, Ostien JT, Salinger AG (2013b) A stabilized assumed deformation gradient finite ele- ment formulation for strongly coupled poromechanical simulations at finite strain. International Journal for Numerical and Analytical Methods in Geomechanics, 37(16):65-75.  

66.Terzaghi Kv, Rendulic L (1934) Die wirksame flachenporositat des betons. Z O ́ st Ing-u ArchitVer 86(1):2 

67.Tezduyar TE, Osawa Y (2000) Finite element stabilization parameters computed from element ma- trices and vectors. Computer Methods in Applied Mechanics and Engineering 190(3):411–430 

68.Truty A, Zimmermann T (2006) Stabilized mixed finite element formulations for materially nonlin- ear partially saturated two-phase media. Computer methods in applied mechanics and engineering 195(13):1517–1546 

69.Wan J (2002) Stabilized finite element methods for coupled geomechanics and multiphase flow. PhD thesis, University of Texas at Austin. 

70.White JA, Borja RI (2008) Stabilized low-order finite elements for coupled solid-deformation/fluid- diffusion and their application to fault zone transients. Computer Methods in Applied Mechanics and Engineering 197(49):4353–4366 

71.Wriggers P, Reese S (1996) A note on enhanced strain methods for large deformations. Computer Methods in Applied Mechanics and Engineering 135(3):201–209 

72.Zhang H, Sanavia L, Schrefler B (1999) An interal length scale in dynamic strain localization of multiphase porous media. Mechanics of Cohesive-frictional Materials 4(5):443–460 

73.Zienkiewicz OC, Chan A, Pastor M, Schrefler B, Shiomi T (1999) Computational geomechanics with special reference to earthquake engineering. Wiley Chichester, UK 









Hi Steve,

Your discussion is interesting. I have included a recent paper that I have found through a web search:

Y. Xie+, G. Wang*. “A Stabilized Iterative Scheme for Coupled Hydro-mechanical Systems Using Reproducing Kernel Particle Method”, International Journal for Numerical Methods in Engineering, Vol. 99 (11), 819-843. DOI: 10.1002/nme.4704.

A link to the paper can be found here:

WaiChing Sun's picture

I also want to add the FEM-FVM coupling model derived by Pr. Prevost's research group at Prinecton and some of its application on CO2 storage problem. This type of multi-model coupling approach is useful for the following reason --- many commerical or research reservoir simulators are written via a finite volume approach, while solid solvers are often written in finite element. While one can always "reinvest the wheel" to make a pure FVM or FEM multiphysics solver, it might be easier (and cheaper) to establish a FVM-FEM coupling if the uncoupled solid and fluid solvers are available to the modelers.   

“Two-way coupling in reservoir–geomechanical models: vertex-centered Galerkin geomechanical model cell-centered and vertex-centered finite volume reservoir models”. J.H. Prévost. International Journal for Numerical Methods in Engineering. 2014. DOI: 10.1002/nme.4657.

The second paper is on studying the thermal stresses of the caprock. Notice that when CO2 is injected underground, it is often injected in its super-critical form, i.e., a combination of pressure and temperture that keeps the CO2 in s dysyr midway between a gas and a liquid. Of course, the CO2 will not have the same temperture of the acquifer or the cap rock and that temperature difference and the bulid up of pressure plume due to the injection may be sufficient enough to affect the integrity of the cap rock.

One thing I notice in this paper is that there is a structural heating term in the balance of energy from both the solid and fluid constituent, which is neglected by numerous previous work. This of course will make the coupling easier to model as it reduces the THM from a three-way coupling problem to a two-way coupling one (i.e. one can solve for the temperature first and then update the displacement and pore pressure without losing accuracy). I also wonder what is the advantage of putting the degree of saturation as the nodal solution? It seems like it may make the THM problem more complicated especially considering the fact that here may be a "degree of saturation jump" between layrs that have profoundly different wettiability. 

“Effect of thermal stresses on caprock integrity during CO2 storage”. G.Y. Gor, T.R. Elliot, J.H. Prevost. International Journal of Greenhouse Gas Control. 12. 2013. pp. 300-309. DOI:10.1016/j.ijggc.2012.11.020.


Jinxiong Zhou's picture

Dear WaiChing, Thank you very much for posting this very interesting topic. A related topic is the modeling of swelling of polymer gels, where the theory of poroelasticity is also used. For gels, we have, typically, two sets of indepedent variables at each node, the chemical potential of fluid and the displacement. I think most of the scheme and procedure for THM is also applicable for modeling gels, albeit some modification is needed. I would like to have your comment on this point.

BTW, can you explain in basic language what is inf-sub condition? What's its effect on numerical stability?


WaiChing Sun's picture

Hi Jinxing, 

Yes, I do think that a lots of the numerical schemes developed for THM problem are appliable to other multiphysics problems. Of course, the subtle difference for various problems are important. In fact if we look at the drained and undrained split slovers by Armero, it really resembles the adiabatic and isothermal solvers developed by Simo and Miehe in the 90s. 

Now regarding the inf-sup condition (often referred as Ladyzhenskaya-Babuska-Brezzi or LBB condition), it is basically a necessary condition for spatial stability. Let me try to explain it in the most simplistic way I could. Note that there are many papers discussing this topic --- I do not think I am the most qualified one to answer this question and hence I would appreciate it if any audience/reader can point out any error or mistake in my explanation. 

First, consider a mutliphysics governing equation that leads to a matrix form that reads (sorry I dont know how to write the xml):

[A B^T ; B C]   [u; p] = [f g]

For poromechanics probelm, the most difficult part is when the pore fluid diffusion is negligible (i.e. near the undrained limit), in which case, we have (assueming incompressible constituents),

[A B^{T; B O] [u; p] = [f 0]

In this case, the material is undrained and the pore-fluid constituent will provide an incompressible constraint to the solid skeleton. One thing we need to pay attention is that B matrix that controls the two-way coupling is not a square matrix but a rectangular one. Now, assume that A is positive definite. We can use the Schur complement B A^{-1} B^T  to solve for p, i.e. ,

B A ^-1 B^T p = B A^-1 f

However, if the null space of (B^T) is spanned by any nontrivial basis,  then there may exist a column vector with non-zero components (let's called it lambda) that satisfies A^-1 B^T lambda = 0. That may cause spatial stability issue because a  spatial oscillatiing (checkerboard pattern) pore pressure field may occur (due to the existence of such a column vector), but we have no "control" over it (i.e. u = A^-1 B^T lambda is zero but lambda is not). In other words, the entire block system becomes ill-posed. 

It is proved that such oscillation occurs when the pore pressure and displacement are discretized by the same set of basis function for the saturated poromechanics problem or other mathematically similar systems (Darcy's flow with V-p pair, incompressible elasticity with u-p pair, Stokes' problem with v-p pair...etc). There are multiple ways to avoid that. For instance, one can modifies the B matrix or add something in the C matrix to elimate the osillation mode.  

Below is an example in which spatial oscillation occurs near the undrained limit in a poro-elastic beam (RIGHT), and the stable solution obtained via a project based stabilized finite element model (LEFT). I have included more examples in my own papers [62,63] and there are many more in the literature. 

I must confess that my explanation is quite shallow. Below is a sample of work concerning the inf-sup condition:

Bathe, Klaus-Jürgen. "The inf–sup condition and its evaluation for mixed finite element methods." Computers & structures 79.2 (2001): 243-252.

Dohrmann, Clark R., and Pavel B. Bochev. "A stabilized finite element method for the Stokes problem based on polynomial pressure projections." International Journal for Numerical Methods in Fluids 46.2 (2004): 183-201.

Fortin, Michel, and F. Brezzi. Mixed and hybrid finite element methods. New York: Springer-Verlag, 1991. 








Jinxiong Zhou's picture

Dear WaiChing, Thanks for your detailed explanation! I learn a lot through reading your post. Jinxiong

WaiChing Sun's picture

Hi Jinxiong, 

I want to add one more point. While inf-sup conditon implies that the matrix B^{T} is injective, injective of B^{T} does not imply that inf-condition holds. In fact, when we consider stability of a given numerical method, we need to consider a sequence of matrices A and B corresponding to a sequence of meshs with mesh size approaching zero. The inf-sup condition is a necessary (but not sufficient) condition that requires the exisits of a positive constant beta that is independent of meshsize h, such that 

x^{T} B^{T} y \geq beta ||x|| ||y||

This is a more strict requirement than the injectivity of B^{T} 





WaiChing Sun's picture

Another new application for inf-sup condition is for concurrently coupling models across different scales. In particular, Bauman et al 2008 has conducted a detailed analysis on the inf-sup condition of  Arlequin method that couples particle and continuum models in 1D. 

Bauman, P. T., Dhia, H. B., Elkhodja, N., Oden, J. T., & Prudhomme, S. (2008). On the application of the Arlequin method to the coupling of particle and continuum models. Computational mechanics42(4), 511-530.

In addition to the mathematical analysis, the inf-sup condition can also be checked numerically via the inf-sup test. The inf-sup test is just a generalized eigenvalue problem corresponding to the inf-sup condition. I have done some work with Dr. Alejandro Mota from Sandia in which we construct an inf-sup test for a multiscale model that couples local and non-local constitutive responses. 

Sun, W., & Mota, A. (2014).  A multiscale overlapped coupling formulation for large-deformation strain localization. Computational Mechanics,54(3), 891-891.

It is interesting to note that both multiphysics AND multiscale problems might lead to block system that requires the finite dimensional spaces used to interpolate different variables being compatible and thus the inf-sup condition.  



Dear WaiChing, Thank you for this very interesting theme. I wanted to let you know that a paper reviewing recent advances in THM fracture will be appearing soon in Advances in Applied Mechanics. Best regards,


WaiChing Sun's picture

Hi Stéphane,

Thanks for pointing out the article to me. I would love to take a look when the book is available.

BTW, I wonder if you have encounter similar inf-sup condition issue in isogeometric models especiually when fomrulating something that leads to a block system (e.g. incompressible elasticity with displacement and hydroastatic pressure as nodal solution)? Of course the basis functions are different between the conventional FEM and isogeometric Galerkin. Any important literature we should be aware of?

Best Regards,



Subscribe to Comments for "Journal Club Theme of September 2014: Numerical modeling of thermo-hydro-mechanical coupling processes in porous media"

Recent comments

More comments


Subscribe to Syndicate