Non-isothermal two-phase flow in deformable porous media: systematic open-source implementation and verification procedure

We propose a formulation for non-isothermal two-component two-phase flow through deformable porous media. The approach covers phase transitions among both phases, i.e. liquid phase components evaporate into the gas phase while gas phase components dissolve or condensate into the liquid phase. These phase transitions always take place in thermodynamic equilibrium. The set of model equations is thereby largely independent of the specific constitutive relations. Starting from general equilibrium equations, we show the evolution of the system of weak formulations of all governing equations, which are then discretised with Taylor-Hood elements in a standard finite element approach. The model equations and the construction of the constitutive equilibria are implemented in the open-source simulator OpenGeoSys, which can be freely used and modified. To verify the implementation, we have selected a number of complementary test cases covering a wide range of process couplings. The numerical model is compared with analytical and semi-analytical solutions of these problems as well as with experimental results. It is shown in the paper that by including thermodynamic effects, solid mechanics, and phase transition processes, the proposed numerical model covers many characteristic features of unsaturated geomaterials and can be employed for the description of a broad range of problems encountered in geotechnical engineering. An open-source FEM tool for non-isothermal two-phase flow in deformable porous or fractured media is presented in detail. The model features phase transitions across both fluid phases based on simple equilibrium conditions. A variety of benchmark tests is presented and compared to other software results and to exact solutions. An open-source FEM tool for non-isothermal two-phase flow in deformable porous or fractured media is presented in detail. The model features phase transitions across both fluid phases based on simple equilibrium conditions. A variety of benchmark tests is presented and compared to other software results and to exact solutions.


Introduction
Numerical models of multiphase flow have been used since the middle of the last century in reservoir engineering and environmental or hydrogeological research, where problems essentially focused on water and non-aqueous liquid interactions. Due to the complexity of the physical processes, those early models were often developed from empirical equations and observations (Hassanizadeh and Gray 1979a). With the upcoming theory of mixtures in the 1970s (Bear 1972;Atkin and Craine 1976;Bowen 1976), multiphase formulations were more frequently based on theoretical approaches, such as continuous fields methods or phase-averaging methods, that describe phases as overlapping continua occupying the entire control volume (Hassanizadeh and Gray 1979b). Today, most numerical multiphase models are founded on continuum mechanical principles (Ehlers 2002;Ehlers and Blome 2003;De Boer 2012), although discontinuous approaches or hybrids of both methods do exist (Pan et al. 2016).
The main principle behind multiphase methods in the framework of porous or fractured media is the consideration of simultaneous flow of at least two immiscible fluids such as water, air, or variations of non-aqueous liquids. The continuum approach purposely ignores the exact inner boundaries between fluid and solid phases and instead assigns volume fractions to describe their continuous distribution over a defined control volume.
Especially if a gaseous phase is to be considered, thermodynamic effects must be included in the model development, since compression and expansion of gases are thermodynamic processes that are often accompanied by temperature changes. Moreover, basic thermodynamic considerations are crucial for capturing phase transition effects, which have a significant impact on pressure predictions in deforming porous media whose pores are filled by both liquids and gases. If components of the fluid phases undergo a phase change due to condensation, dissolution or evaporation, the corresponding energy cost must be considered in the energy balance of the multiphase body.
Furthermore, widely different fields of engineering tasks require the consideration of solid mechanics in order to evaluate structural safety, consolidation phenomena or the change of the porous solid's properties in regard to fluid transport. This requires to account for the rate of change of linear momentum to be included in the set of balance equations and for an appropriate constitutive setting based on the generalised Biot theory (Biot 1941).
Mechanical models that account for fluid flow and heat transport are often referred to as thermo-hydromechanical (THM) models. Several THM codes are available in the literature. Many of them refer only to fully saturated conditions (Xu et al. 1999;Ai and Wang 2016;Zheng et al. 2017;Cui et al. 2018), to unsaturated media with an isobaric gas phase (Sheng et al. 2003;Rutqvist et al. 2008;Wang et al. 2009Wang et al. , 2015Zhu et al. 2020), or to formulations that do not cover phase transitions (Schrefler and Xiaoyong 1993;Oettl et al. 2004;Klubertanz et al. 2003).
There is also a variety of models available that exclude mechanical processes (Pruess and Joseph 1984;Class et al. 2002;Pruess 2004;Bourgeat et al. 2013), or those that consider isothermal conditions only (Khoei and Mohammadnejad 2011;Dagher et al. 2019).
The isobaric gas phase approach uses the Richards approximation 1 . In this approximation, the gas phase pressure is considered constant and bulk gas flow is assumed to be instantaneous and thus ignored so that no mass balance equation for the gas phase is required. This approximation produces reasonable results only when gas permeability is large enough and the gas phase can drain freely (Wang et al. 2011). This statement indicates that the isobaric-gas approach is not suitable for tasks involving the buildup of gas pressures or the evaluation of barrier performance when subjected to significant gas pressures. Such phenomena are of interest in several deep geotechnical applications such as geothermal reservoirs or repositories for nuclear waste disposal. For such applications, a flexible and fully tested modelling approach is advantageous, which can assess all phases of deep geological repositories for various geotechnical applications, i.e. their construction, operation and post-operational phases, within the same numerical setting.
Page 3 of 29 107 Vol.: (0123456789) Non-isothermal multiphase flow models, which take phase transitions and the deformation of the solid phase into account, are not often encountered so far. Examples of such simulators are implemented in the model environments CODE_BRIGHT (Olivella et al. 1996), MOOSE (Green et al. 2021) or are available as a combination of the simulator TOUGH and the commercial code FLAC 3D (Rutqvist 2011). Although these examples provide similar functionality to the model presented here, they have some disadvantages, either because they are not freely usable or because they do not have open source code, offer only a limited choice of constitutive relations, or because they only provide sequential couplings between processes, which are partly based on the exchange of text files.
In this paper, we propose a TH 2 M-model, where H 2 indicates the two-phase character of our THM formulation. The non-isothermal formulation is founded on the theory of porous media, prominently developed by De Boer (2012); Ehlers (2002), and considers a triphasic medium consisting of a deformable and intrinsically compressible solid skeleton whose pores are filled by two immiscible fluid phases. Phase change and phase transitions from one fluid phase to another are permitted via evaporation, condensation, as well as dissolution of fluid constituents. Thus, this model approach is a two-phase, two-component formulation. Both components can be transported in phase space, but are subject to the constraints of phase equilibrium.
Foremost fields of use for this model setup are geotechnical applications, such as radioactive waste disposal, carbon dioxide storage in geological formations, enhanced oil or gas recovery, or deep geothermal energy systems. In the context of nuclear waste management, the model is currently utilised for benchmark studies and the design of underground rock laboratory experiments (c.f. the DECOVALEX project 2 and the EURAD program, 3 BenVaSim project, cf. Lux and Rutenberg (2018)).
The consideration of solid mechanical processes is required to assess barrier integrity, soil consolidation or changes in host rock porosity or permeability in underground storage applications. It also allows the inclusion of shrinking and swelling processes and their effects on stress and displacement fields.
The inclusion of phase transition processes provides knowledge of phase compositions, which permits detailed investigations of physical phenomena such as vapour diffusion or buoyancy-driven advective flow. For example, the dissolution of gases (e.g. hydrogen) in porewater is one of the main processes responsible for the reduction of gas pressure in radioactive waste repositories. Phase transition processes are also important for geomechanical analyses, as otherwise fluid pore pressures (and thus effective stresses) may be quantified incorrectly (Khaledi et al. 2021). It must be noted here that our formulation does not cover the appearance or disappearance of individual phases. Such tasks require specific methods such as the switching of or the definition of persistent primary variables (Lauser et al. 2011;Bourgeat et al. 2009;Mahjoub et al. 2018).
In addition to increasing functionality of codes for numerical simulation of coupled multi-field processes in fractured and porous media, the availability of mature open-source platforms as a starting point for further development and application is becoming a community effort. Open-source solutions are important for both flexibility and transparency. Flexibility is a prerequisite for workflow integration, i.e. seamless connection of simulation tools into data integration and data analytics frameworks . Providing transparent simulation tools is fundamental when it comes to improving public acceptance of large-scale geotechnical projects which are often controversially discussed. Hence, since siting decisions are usually made by politicians on the basis of scientific suggestions, transparency is crucial to increase public confidence in scientific methods. Open-source workflows require significant extra effort for quality assurance on the side of the scientific community concerning software engineering including version management, automated benchmarking, peer-review for code development etc. The presented TH 2 M model is implemented into the OpenGeoSys (Bilke et al. 2019) framework and permanently available through the related GitLab repository. 4 OpenGeoSys offers a variety of coupled processes, including those featuring multiple fluid phases. However, the presented paper offers the first monolithically coupled model of nonisothermal two-phase flow with mechanics The structure of this article is as follows. As part of the introduction, we summarise all assumptions made during the development of the formulation.
In Sect. 3, the mathematical framework based on continuum mechanical balance equations is outlined, followed by a detailed description of the constitutive setting for phase equilibria as well as thermodynamical, mechanical and transport properties in Sect. 4. In Sect. 5, the balance equations are developed further to governing equations by applying further assumptions and by selecting primary variables. Then, weak formulations of the governing equations are obtained by applying a standard Bubnov-Galerkin procedure in Sect. 6. Section 7 describes briefly some aspects of the numerical implementation, while Sect. 8 presents a selection of benchmark tests for code verification purposes. Finally, Sect. 9 concludes and summarises the article and suggests necessary future work.

Notation and assumptions
The work at hand describes a triphasic body consisting of a deformable and intrinsically compressible solid skeleton phase = S as well as two compressible, interstitial fluid phases ∈ {L, G} . The fluid phases are composed of the constituents ∈ {W, C} that are allowed to cross the phase boundaries according to specified phase-equilibrium conditions. To conveniently describe the phase transitions, the fluidrelated equations are developed in terms of these constituents rather than actual phases. Therefore, the fluid constituent equations are perfectly symmetrical, since both constituents are allowed to cross phase boundaries and both are considered compressible.
We chose gas phase pressure p GR and capillary pressure p cap as primary variables for the hydraulic part of the formulation. While pressure-saturation approaches are often used since they allow for an easy assignment of boundary conditions, they are known for their inability to reproduce physical behaviour in heterogeneous media (unless upwind correction techniques are employed). Pressure-pressure formulations allow for correct physical behaviour in this context without such correction (Park et al. 2011). For mechanical and thermodynamical processes, we use displacement of the solid phase u S and temperature T as primary variables.
Multiphase fluid flow, solid deformation and temperature evolution affect each other strongly. Due to its high accuracy and numerical stability, we chose a monolithic coupling concept where all coefficients of the partial differential equations are assembled into one global matrix-vector system and solved simultaneously.
Further assumptions that form the basis of our formulation include: -solid deformation follows the small-strain approach, where absolutes of the displacement gradient remain below about 5 % , so that linear kinematics can be used; -phase motion occurs in a quasi-static manner, a = 0. -local thermal equilibrium, T = T is maintained; -all phases are compressible according to suitable equations of state; -phase transitions occur among fluid phases only; -inter-phase mass transfer occurs instantaneously and is thus subject to equilibrium conditions;

Balance equations
Considering a continuum mechanical approach, a multiphase system can be expressed as a mixture of constituents. According to mixture theory (cf. Green and Naghdi 1969;Passman et al. 1984), balance equations for an entire multiphase body can be assembled in the same manner as for individual constituents, insofar all interactions between all involved constituents are incorporated. Therefore, general balance equations can be taken directly from continuum mechanical methods of single phase bodies (Ehlers 2002). A generic balance equation relates the temporal change of a quantity within some reference volume with effluxes over its boundary, supply resulting from external distance, and production due to interaction with other phases or constituents (internal interactions of the mixture). Which mechanical or thermodynamical quantities are necessary to be considered depends on the purpose of the model and its assumptions. In this work, it is sufficient to incorporate the Page 5 of 29 107 Vol.: (0123456789) balances of mass, linear and angular 5 momentum as well as the energy of the multiphase body.

Mass balance
The local form of the mass balance for an arbitrary phase can be written as where (•) � denotes the material time derivative of a quantity • with respect to the phase velocity of , and ̂ are density production terms related to mass transfer from or to phase . Here, we consider only physical inter-phase mass transfer of the constituents, i.e. the latter undergo no chemical alteration. Thus, the phase density production is given as with as phases other than and with as constituents residing in all phases. Indices define the direction of mass transfer from phase to phase .

Solid phase
The solid phase consists of a single component S only, while mass exchange to other phases is precluded by assumption, such that and for simplicity, it is convenient to write S S = S . Since there is no mass production, the mass balance equation for the solid phase simplifies to where u S as the displacement vector of the solid phase.

Component mass balance equation
The mass transfer terms ̂ in (1) express the combination of all phase transitions of all constituents. In case that consists of multiple constituents (as in this work), it may be convenient to devise balance equations based on constituents within , which yield with component velocity v = v + d , which results from phase velocity v and diffusion velocity d . Except for the solid material, all components considered in this work can exist in both fluid phases ∈ {L, G} . To describe a component balance for the entire system, the mass balance equations for a balanced component are summed over all fluid phases. Since all phase transitions are restricted to fluid phases, the mass transfer terms disappear since mass fluxes of leaving one phase are equal and opposite to those entering the other phases.
When all material time derivatives are transformed into a solid phase reference frame, the mass balance for a constituent reads with the diffusive-dispersive mass flux J = d and the advective mass flux where w S is the phase velocity relative to the movement of the solid phase. Continuing evaluation of the balance equation is aided by the introduction of phase saturation s and phase volume fraction = s , where phase saturation is the ratio of phase volume to the combined volume of all fluid phases and is porosity. Phase saturation and porosity, together with the real (intrinsic) mass density, define the apparent mass density = s R . Temporal changes of are caused by the combined changes of all of those quantities.
The solid phase mass balance (4) can be inserted into (6) in order to eliminate porosity variations, which finally leads to Note in passing, that this commonly used form of the mass balance already includes the assumption of a variable porosity-as can be observed by the different factor in front of the solid velocity divergence between Eq. (6) and Eq. 7-and requires specification of an equation of state for the intrinsic solid density evolution.

Balance of linear momentum
The variation of linear momentum of the multiphase body is caused by the combination of all acting forces. Thereby, it can be differentiated between body forces acting on the control volume (e.g. gravity) and forces acting on its surface. The local form of the balance of linear momentum can be written as The partial stress tensor results from the partial surface traction vector t following Cauchy's fundamental theorem t = n . A momentum production term ŝ has to be considered for a phase-formulation of the momentum balance equation, which results from direct momentum productions p and from material phase transition momenta ̂v . The production terms follow the constraint Considering this constraint, the final balance equation of linear momentum can be written as where a = v � is phase acceleration and b are specific body forces.

Energy balance
The total specific energy per unit mass is composed of the specific internal energy and the specific kinetic energy 1 2 v ⋅ v . The rate of change of total energy implies a change in both internal and kinetic energy and is equal to the work done by external body forces and surface tractions, the heat input by conduction and external supply, and energy production terms (e.g. due to phase transitions).
The local form of the total energy balance can be written as follows with specific internal energy u . Since only small velocities are expected, all parts associated with kinetic energy are neglected. The term governing the power of the surface tractions can be expanded to with the spatial velocity gradient l = grad v and the deviatoric part of the rate of deformation tensor d D . By incorporating the momentum balance equation (10) and neglecting deviatoric stress power contributions, we find Rewriting time derivatives to follow the motion of the solid phase and introducing specific enthalpy 6 h = u + p −1 , the final balance equation reads In analogy to Sect. 3.1.2, the above equation can also be set up for arbitrary components in the phase . One obtains with the relative component velocity w S = v − v S considered relative to the solid phase.

Constitutive setting
The system of balance equations from the previous chapter has to be closed by a variety of constitutive relations, which in turn specify the necessary material properties and express the thermodynamic equilibrium between the constituents of liquid and gas phases. The specific choices of constitutive equations made in this section are used in the examples shown in Sect. 8, but the process equations are derived independently so that a large variety of material laws and equations of state can be utilised.
In order to describe the phase composition, we assume both fluid phases to consist of a binary mixture of components W and C , where W is considered the main component of the liquid phase L (e.g. water) and C is considered the component mainly present in the carrier gas phase G.

Fluid phase composition
As stated above, both constituents W and C exist in both fluid phases; their distribution among the phases follows a very basic equilibrium condition which is based on the ideal gas conception and Henry's law. A more sophisticated phase equilibrium model for twophase flow applications is proposed in Grunwald et al. (2020).

Gas phase
Determining the composition of the gas phase, we assume that the sum of all constituents' partial pressures accounts for the entire gas phase pressure (Dalton's law). Therefore, the binary gas phase composition in terms of molar fractions x n,G results in with the partial pressure of water vapour p W GR , which is derived from the saturation vapour pressure of pure water p W vap . The saturation vapour pressure can be determined by the approximate Clausius-Clapeyron equation, Vapour pressures derived from above-mentioned methods refer to flat liquid-vapour interfaces only. Above curved interfaces such as those prevailing in porous media, the vapour pressure may change depending on the direction of the curvature. In general, convex interfaces as on droplets have larger vapour pressures, concave surfaces exhibit smaller vapour pressures than flat interfaces.
The curvature of the vapour-liquid interface in porous media is directly associated with wettability and the capillary pressure. Written in terms of capillary pressure, the Kelvin-Laplace equation can be used to evaluate the true vapour pressure above curved interfaces: with M W and W LR as molar mass of component W and real mass density of the liquid phase. (19) shows that this effect is significant for very large capillary pressures only (to lower the vapour pressure by 1 % requires a capillary pressure of approx. 1.3 MPa ), as they are characteristic for media with very small pores such as claystone.

Liquid phase
The concentration of component C dissolved in the liquid phase is proportional to its partial pressure in the gas phase with Henry-coefficient H C . A relation for Henry coefficients as a function of temperature can be found in Sander (2015), which reads where H is given in mol m −3 Pa −1 .

Gas phase
The gas phase density can be determined using an appropriate equation of state for binary mixtures, such as the Peng-Robinson equation of state (Peng and Robinson 1976). For simplicity, we use the thermal equation of state of ideal gases in this work, which gives sufficient results at high temperatures and low pressures. A more realistic approach based on realgas equations of state is developed and described in detail in Grunwald et al. (2020). The ideal gas law for the entire gas phase and for gas phase constituents can be written as with the average molar mass M G = ∑ x n,G M and constituent mass fraction x m,G = GR −1 GR .

Liquid density
As a model for the liquid phase density, a simple, empirical multi-linear equation of state can be fitted easily to measurement data in a state region of interest: where c,LR is a liquid phase coefficient to describe density changes due to concentration changes of component C . With LR = Σ =W,C LR and W LR = LR p LR , T, c C L = 0 , we find the partial density of dissolved gas in the liquid phase via (20) and (23) to

Specific enthalpy and specific internal energy
Specific enthalpies of the fluid phases are considered as functions of composition, realised as the mass fraction weighted sum of the constituents' enthalpies, with Δh W L = Δh C G = 0 , Δh W G = Δh W vap is enthalpy of vaporisation, and Δh C L = Δh C solv is the enthalpy change of solution.
For determining the internal energy, contributions of the specific phase enthalpies and (for fluids) the pressure work are added: In the case of the solid, this contribution is ignored and u S = h S applies. For a detailed discussion on this topic refer to Lion et al. (2017), for example.

Heat conduction
Heat conduction is driven by the temperature gradients and proportional to the heat conductivity tensor . This relation is known as Fourier's law and yields When heat conduction occurs over multiple phases, a mixing rule can describe averaged heat conduction if local thermal equilibrium is assumed. In this case, we apply a very simple model (upper Wiener bound, Wiener (1912)) to find an effective heat conductivity by averaging individual phase conductivities by volume fraction, which resembles a parallel conductivity model where the overall conductivity is the arithmetic mean of individual phase conductivities. Other thermal conductivity mixing rules are possible and the choice of the model strongly depends on the materials and material distributions present in the model domain. A large variety of proposed mixing models for thermal conductivity of unsaturated porous media is reviewed in detail by Dong et al. (2015) and Tong et al. (2009) Vol.: (0123456789) provides multiple mixing rules for multiphase flow problems in particular.

Molecular diffusion
The diffusive-dispersive mass flux J results directly from the product of the component partial density R and the diffusion velocity of the component d = v − v . This diffusion velocity denotes the difference between component velocity and phase velocity. In addition, Fick's law describes the diffusive mass flux as a function of the concentration gradient. Given in form of mass fraction gradient (Hassanizadeh 1986a) it reads with the effective molecular dispersion tensor D ,m = D . In Hassanizadeh (1986b), the general case of Fick's law is described as so that describes the diffusive-dispersive velocity of a component in the phase . More sophisticated multicomponent gas diffusion models are discussed in e.g. Pisani (2008).

Viscosity
Dynamic viscosity of composite phases is obtained from the simple mixing rule which computes the molar fraction averaged combined viscosity R from constitutent viscosities R . More elaborate, classical mixing models can be found in Wilke (1950) for gases and in Grunberg and Nissan (1949) for liquids.

Saturation
Phase saturation s is the volumetric fraction of phase within the pore space. It is considered in this work as a function of the capillary pressure only. Depending on the field of study, there are numerous different correlations available that correlate liquid phase saturation with capillary pressure, most prominently the works of Brooks and Corey (1964) or Van Genuchten (1980), which both are included in the material property library of our modelling platform. For more complex problems, the use of more sophisticated models is advisable. A detailed review of the existing soil properties can be consulted in Leong and Rahardjo (1997).
The Brooks-Corey formulation accounts for residual saturations of both fluid phases and links effective saturation to capillary pressure by where p b is the entry pressure (i.e. the minimum capillary pressure for desaturation to occur) and is a parameter governing the shape of the curve. This effective saturation is defined by with residual saturation s res L and maximum saturation s max L = 1 − s res G . Since only two fluid phases are present, it follows that which constitutes the fundamental saturation constraint.

Relative permeability
Due to the existence of multiple phases within the pore space, the movement of a fluid phase is obstructed by the presence of the other phase. In multiphase flow applications, this effect is usually realised by introducing relative permeabilities as functions of saturation which calculate the effective permeability of each phase as described in the extended Darcy law.
Relative permeability models can be purely empirical or derived from theoretical pore network models. The latter have the advantage that they can be differentiated on the basis of capillary pressuresaturation relationships (Helmig 1997). Several formulations for relative permeability are included in OpenGeoSys' material property library, such as Mualem (1976) or Brooks and Corey (1964) which both are based on the theoretical approach. The Brooks-Corey models for wetting and non-wetting phases are and with empirical constant and entry pressure p b from the Brooks-Corey characteristic curve (33). Once more, the choice of the appropriate relative permeability model depends strongly on the problem under consideration.

Solid phase density and porosity
Assuming solid phase density being a function of temperature and the hydrostatic pressure in the solid, we recognize the latter is composed of pore fluid pressure p FR = s G p GR + s L p LR and an effective pressure p E resulting from effective stress, The solid density rate follows where K SR is the intrinsic bulk modulus of the solid grains and K Sp is the compression modulus of the grains measured under a change in pore pressure in unjacketed conditions. Now we split the total pressure rate into fluid and effective pressure rates. We further assume that volume changes associated with homogeneous temperature changes and homogeneous pore pressure changes (as in the aforementioned unjacketed test) do not cause a change in effective stresses: where K S is the bulk modulus of the solid skeleton. This equation can be rearranged to introduce a modified effective stress that no longer does include an explicit dependency on the pore pressure 7 . Instead this effect is captured by the introduction of Biot's coefficient B : . Total pressure and real solid pressure are equal in the unjacketed test while the effective pressure remains zero. Hence, experimentally we get K Sp = K SR . The second term in (39) is treated by recognising that p SR ∕ p E = (1 − ) −1 and by substituting the definitions in (40): with p,SR = K −1 SR . From solid phase mass balance (4) we find for the rate of porosity change. (40) Page 11 of 29 107 Vol.: (0123456789)

Other relations
The compilation of constitutive equations is not meant to be complete here as this would go beyond the scope of this paper and is available in the cited textbooks above. We focus here on the constitutive theory important for the examples in Sect. 8 and from the viewpoint of code implementation presented in Sect. 7. Some other essential and mostly standard constitutive relations will be given without detailed discussion in Sect. 5.

Governing equations
In this section, balance equations are specialised to the triphasic system ∈ {G, L, S} and the independent unknowns are chosen to be gas pressure p GR , capillary pressure p cap , temperature T, and displacement u S .

Mass balance
From (7), the constituent mass balances are added for each fluid phase ∈ {G, L} , where the solid phase density derivative is replaced by (42). Expanding the remaining pore fluid pressure derivative leads to with the abbreviations A = A L + A G for combined constituent advective fluxes, J = J L + J G for combined constituent diffusive-dispersive fluxes, and FR = s G GR + s L LR as an intrinsic effective constituent mass density. The dispersive-diffusive mass fluxes are driven by gradients in constituent mass fraction x m, . In this work, the distribution of a constituent among phases is governed by a phase equilibrium model fully determined locally by the thermodynamic state of the system. Therefore, the mass fractions are understood as functions of the independent unknowns that define the state of the multiphase body and the diffusive-dispersive mass flux can be expanded as with binary diffusivity defined in Sect. 4.5.

Momentum balance
The phase stress in the general momentum balance equation (10) can be split into parts contributed by the pore pressures and extra stresses with E ≈ 0 for fluid phases ∈ {L, G} due to the insignificance of the associated forces in comparison to pressure gradients and frictional interaction with the solid (Markert 2006). Written in terms of fluid phases only, (10) reads into which an extra momentum production term In (48), k rel s is a function of phase saturation and k S is the intrinsic permeability of the porous medium.
Introducing Biot's coefficient to account for the solid grain compressibility, we follow the general approach in Lewis and Schrefler (1998) to find the (45) for ∈ {L, G} 8 Note that an extra term p cap s −1 L grad s L appears in the relation for the liquid's Darcy velocity. This term is neglected here, see (Häberle 2017) and references therein for a discussion. effective stress (see also the discussion in Sect. 4.9) for the solid phase to where p FR = s G p GR + s L p LR is the overall pore fluid pressure. We consider the linear momentum balance of the entire multiphase body by summing up the balance equations for all phase momenta. From the constraints in (2) and (9) as well as assuming b ≡ b , the overall aggregate momentum equation reduces to with an overall effective density = S SR + L LR + G GR .

Energy balance
Following the assumption of local thermal equilibrium, the overall aggregate energy balance is obtained by summing (14) for all phases and components which cancels out all energy productions ( ∑û = 0 ). Incorporating the momentum balances for fluids (47) we find p = grad p − b and by considering (9) and (25) we obtain the overall aggregate energy balance with with the combined effective heat conduction q = ∑ q = − ef f grad T governed by Fourier's law as described in Sect. 4.4.

Weak formulations
In this section, the balance equations for mass, linear momentum and energy are restated into weak forms by multiplication with a Bubnov-Galerkin-type test function a and integrated over the entire domain Ω.

Component mass balance
In this work, two constituents are considered where each component can exist in both fluid phases simultaneously. Therefore, (44) is applied for each constituent individually. The mass flux terms can be expanded using partial integration   (16) and Henry-coefficient in (20) to zero), (54) is reduced consequently to the classical mass balances of gas and liquid phase.

Momentum balance
Similar to the previous section, the linear momentum balance equation (50) is now integrated over the overall volume Ω and multiplied by a vector-valued test function u s . Substituting the total stress tensor defined in Sect. 5.2 and applying the divergence theorem, we find the weak formulation for the momentum balance to with the total traction vector t and with s = sym grad u s .

Energy balance
Similar to mass and momentum balance equations, the energy governing equation, (51), is integrated over the domain Ω and multiplied by a test function T . Conductive heat fluxes are governed by Fourier's law (28). The advective and conductive-diffusive heat fluxes split into internal heat transport and heat transport across the domain boundary Γ . By applying the divergence theorem we find

Implementation
The weak forms of constituent mass balance, (54), overall aggregate momentum, (55), and energy balance, (56), are discretized using a standard finite element procedure. Thereby, the field of displacement unknowns u S can be discretized by elements of higher order than the unknowns gas pressure p GR , capillary pressure p cap , and temperature T, utilising Taylor-Hood-elements. For the temporal discretization using a Backward-Euler scheme, a variable time step length Δt is employed. Scalar-valued unknowns a and test functions a , their rates and gradients are approximated via where N a ( ) are row vectors of appropriate shape functions and â are column vectors of nodal degrees of freedom, represented by The entries N i a ( ) are the interpolation weights associated with node i of the finite element at position 9 .
For the description of displacement processes it is advantageous to express vector and tensor values in the form of column matrices while maintaining their (56) tensor character (Nagel et al. 2016). The vector-valued displacement u S and the scalar volume strain rate div u S � S are approximated by with m as the vectorial counterpart of the identity tensor I given by m T = [1 1 1 0 0 0] and with quadratic shape functions N u arranged such that the resulting matrix accommodates for vectorial nodal unknowns. B u is the corresponding strain-displacement matrix; cf. Nagel et al. (2016) for implementational details. Instead of analytically expanding the time derivatives of secondary variables in terms of the independent unknowns, they are approximated here by backward difference quotients directly. Due to the temporal discretization of the numerical method, this has proven to represent more accurate rates than using derivatives only exact in the limit Δt → 0 (Celia et al. 1990). We use the notation with y t+Δt denoting a function evaluation at the current time step, y t at the previous time step and Δt as the time step size.
The discretised forms of the component mass balances, momentum and energy balances are listed in Sect. 1. The resulting system of equations is solved in a monolithic way using an incremental iterative Newton-Raphson scheme.
The set of equations together with the constitutive setting is implemented into the open-source software OpenGeoSys (Kolditz et al. 2012;Bilke et al. 2019), which can be downloaded from the project website 10 . Together with the source code, the website provides an automated documentation of the code and a large variety of test cases and benchmark tests that verify the implementation and can be used as starting points for new users in order to develop own numerical models. The quality of the software is maintained by using state-of-the art development tools such as git, and new contributions to the source code are subjected to a peer-review process.

Verification
Approximate numerical solutions of coupled, nonlinear problems are usually verified by comparison to closed-form solutions. Those exact solutions exist only for a very limited range of spatial domains or for simplified constitutive settings . To our knowledge, there is no closed-form solution available that covers all aspects of non-isothermal two-phase flow in deformable media as described in this work. There are a few examples with roughly the same functionality (cf. LaForce et al. 2014a, b), but these do not consider phase transitions or diffusion processes. Instead, the task of verifying the developed set of equations was split up into several parts. Thereby, the complexity of the equation system is reduced to simpler combinations of the individual physical processes, cf. Fig. 2. Each of the meaningful subsets of process combinations can then be tested against closed-form solutions or the results of physical experiments. An entire set of test cases (benchmark tests) is produced to verify the individual terms in the system of equations presented in previous chapters.
All of these test cases are available on the project website 11 , where they can be downloaded together with the source code of the entire modelling platform. Examples of the first and second levels of simplification, i.e. TH (Sect. 8.1), THM (Sect. 8.2), TH 2 (Sect. 8.3), and H 2 M benchmarks (Sect. 8.4), respectively, are shown and explained in detail in the following sections.

TH: Isentropic gas compression
In this test, the thermodynamic relationships between gas pressure, temperature and density are tested. For that, a cube-shaped domain consisting of an ideal gas is compressed by 50 % of its initial volume over a short period of time, starting from the top surface. The boundaries of that domain are impermeable to fluid flow, therefore gas pressure and density must increase as a result of the decreasing volume. Since energy flow across the boundaries is also prevented, this compression is an adiabatic change of state. No frictional losses are taken into account, so the process Page 15 of 29 107 Vol.: (0123456789) can be reversed at any time and the entropy in the overall system remains constant. The mass balance for such a system can be found by simplifying (7). With = s G = B = 1 and A = J = 0 one obtains or where (e) � S = div u S � S is the volume strain rate. This function can now be integrated over the density interval 0 GR , GR and the volume strain interval [0, e] so that the final evolution equation of the gas density is given by with initial gas density 0 GR . The evolution of gas pressure can be found in a similar way from the energy balance (51). Without fluid velocities, heat conduction, or gravity, the energy balance for the gas in this case simplifies to with u G = h G − p GR with initial temperature T 0 . The numerical approach for this problem is carried out with a very simple setup. The gas cube is represented by a single quadrilateral element, which eliminates any transient fluid velocities or heat conduction effects. The liquid phase equation was disabled by defining s G = 1 as a constant parameter, which also prevents any phase transitions. Due to the simplification of the model, only the specific heat capacity and the molar mass of the gas phase remain as relevant material properties, which are given here as c pG = 14.3 kJ kg −1 K −1 and M G = 0.002016 kg mol −1 and which correspond to molecular hydrogen. Other material properties do not affect the model result, except for those kept constant for manipulation and simplification of the complex TH 2 M equation system: The saturation of the liquid phase remains at s L = 0 , the porosity of the medium and the Biot coefficient of the solid phase are given by = B = 1.
As can be seen from the equations (63), (66), and (67), the analytical solutions are scaled via the initial conditions. Thus, initial pressure and temperature are irrelevant in the numerical model as long as the gas density is consistently determined from the starting conditions via the ideal gas law. Compression is carried out by a displacement boundary condition of u z = −0.5 m at the upper boundary of the domain, which is applied via a time ramp over a period of 10 s . This time period is divided into 100 discrete time steps of Δt = 0.1 s each.
From Figure 3 it can be seen that the relative density when the cube is fully compressed (at t = 10 s ) is about 65 % above its initial value. This fact may seem counterintuitive, after all, the height of the cube is shortened by 50 % . Since no mass can escape, one could expect a doubling of the gas density. Here the effect of the assumption of small deformations made at the beginning of this paper becomes apparent: In this approach, the integration volume remains fixed at its initial value (reference configuration). Geometrical changes are only included to first order (linear kinematics). It can be deduced from the results of this example in which the height of the domain is shifted by 50 % , this basic assumption is no longer valid. Regardless, the analytical solutions (63), (66), and (67) as well as their numerical approximation are valid for the verification of the thermodynamic model. Both solutions agree very well in all cases, which shows that the constitutive approaches of specific internal energy, specific enthalpy and the contribution of pressure work have been implemented correctly.

THM: Consolidation around a point heat source
Containers with highly radioactive waste generate a lot of heat in the repository, which is released into the surrounding rock. A problem considering thermally induced expansion of liquid and solid phases resulting in thermal strain was first proposed by Booker and Savvidou (1985), who provided an analytical solution for temperature, pressure, displacement and stress evolution in a spherical domain around a point heat source. This work was later revised and clarified by Chaudhry et al. (2019).
A point heat source will increase the temperature within the biphasic domain that will cause expansion of the soil water and the solid skeleton. Since pore water expansion exceeds the expansion of the pore space of the skeleton, the rise of pore water pressure will cause a reduction of geomechanical effective stress. Over time, the pore water flow directed outwards will result in the dissipation of the pore water pressure and thereby the relaxation of the effective stresses.
We use the revised analytical solution (Chaudhry et al. 2019) to verify our numerical formulation. Therefore, our model complexity has to be reduced by one degree of freedom. This is realized by constraining capillary pressure to zero over the entire domain, while the gas phase pressure equation will take the role of the pore water. Since the analytical solution is developed for spherical domains, the numerical model domain is set up as a radial symmetrical quadrant, c.f. Fig. 4.
At boundaries x = 0 m and y = 0 m , displacements in y or x direction are restricted while temperature and pressure are kept constant at the outer boundary at r = 10 m . A constant heat source of q th = 300 J s −1 is applied in the center an r = 0 m . All material properties used in this test case are given in Table 1. In Fig. 5, it can be seen that temperature and pore pressure at a point at r = 0.5 m rise over time as expected. The pressure build-up of the pore water causes an increase of the tensile stress in the solid, whose maximum occurs at approximately t = 3.0 × 10 5 s . Overall, all parameters plotted in Fig. 5 show very good agreement with the analytical solution of the problem. The small deviations seen in the pressure and effective stress results are due to the fact that the analytical model contains expressions with an invariant density, while the numerical model updates density consistent with the specified thermal expansivities and compressibilities of the pore fluid. Thus, the deviation is the result of a simplification made in the derivation of the analytical solution but not in the numerical model.

TH 2 : Heat pipe
Similar to the former example, the test resembles a high-level radioactive waste repository generating heat in a low-permeable medium, which is this time considered rigid but of variable saturation. Driven by the heat source in the centre of the domain, liquid water evaporates and thus increases the vapour pressure in the gas phase. The resulting pressure gradient causes a gas flow directed away from the heat source. At the cooler parts of the domain, vapour condenses again and returns, driven by suction, as liquid water towards the heat source.
If the heat source is strong enough, the regions adjacent to the heat source will dry out entirely, while the distant regions are not affected by a change in saturation. In these two border regions, single-phase conditions prevail (gas phase at the heat source, liquid phase in the distant region, cf. Fig. 6). The region in between is a two-phase region and represents the actual heat pipe zone. Here, the cycles of evaporation, gas flow, condensation and liquid backflow take place. In the stationary case, the mass flows of water vapour and liquid water are equal.
This test case was originally proposed by Udell (1985), who present a semi-analytical solution of the problem. This solution is valid for the steady-state and the heat-pipe region only. Departing from the method proposed by Udell (1985), the semi-analytical solution in this work was extended as presented by Helmig (1997) to include ideal gas behaviour and a mixing rule for gas viscosity (cf. (32)). Although a radially symmetric version has been proposed for this purpose by Doughty and Pruess (1988), the classical, Cartesian formulation is used in this work to verify the numerical model. Schematics of the concept of this test are shown in Fig. 6.
A constant heat flow q th = 100 W m −2 is applied at the right model boundary x R = 1 m . Constant pressure and temperature conditions prevail at the left boundary x L = 0 m of the domain. Consequently, the saturation and the molar fraction of air in the gas phase do not change with time at this point. Over time, a saturation profile forms with a gradient from the heated end to the cold end. In conjunction with the relationship between capillary pressure and saturation, this leads to a pressure gradient in the aqueous phase in the same direction. The pressure gradients of the two phases run in opposite directions, so that finally a stationary circular flow is established. The water vapour flows away from the heat source, taking with it the enthalpy of evaporation, which is released again during condensation. Since the condensed water simultaneously flows back to the heat source due to capillary forces, this enables a very efficient heat transfer. The steam pushes the non-water gas component to the left, where parts of it leave the problem domain during the build-up of the steady-state.
Saturation and relative permeability models are governed by the Brooks-Corey formulations model (c.f. Sects. 4.7, 4.8). The model domain was discretized using a one-dimensional mesh consisting of 200 quadrilateral elements with a width of dx = 0.5 cm each. Initial and boundary conditions are shown in Fig. 6 and material properties are summarised in Table 2. Resulting temperature, saturation, and air molar fraction distributions in the domain are shown in Fig. 7 at the steady-state, which establishes after t ≈ 5 × 10 6 s . All three variables reproduce the results of the analytical solution very accurately. Due to the complexity of the processes involved, a successful implementation of the model equations and the constitutive relationships can therefore be assumed.
The length of the resulting heat pipe is approximately 90 cm In Fig. 7, it can be seen very clearly how the increasing vapour pressure caused by the rising temperature elevates the gas pressure. The resulting pressure gradient transports all air present in the domain to the left edge of the model. From approximately x = 10 cm there is almost no more air in the gas phase.
As expected, the temperature in the region does not rise significantly despite the constant heat source, but remains in the range of the boiling temperature of water. This is due to the fact that almost all supplied energy is used for the ongoing evaporation process (latent heat). The temperature difference between the left and right model boundary is about ΔT = 9 K , from which an apparent thermal conductivity of HP = 11 W m −1 K −1 can be derived (compare the individual thermal conductivities of the three phases present in the model area in Table 2). The heat pipe effect ensures that a considerably larger amount of heat is transported through the domain than would be the case with conventional conductive heat transport.

H 2 M: Liakopoulos experiment
Multiphase flow coupled with consolidation is tested using the well-known experiment of Liakopoulos (1964), which has become one of the most frequently consulted benchmark tests for multiphase flow processes in deformable media (cf. Schrefler and Xiaoyong 1993;Helmig 1997;Lewis and Schrefler 1998;Ehlers et al. 2004;Khoei and Mohammadnejad 2011;Asadi and Ataie-Ashtiani 2015). The physical experiment consisted of a narrow acrylic column with a length of 1 m packed with Del Monte sand. In order to guarantee uniform flow conditions, a continuous water supply was added to the top of the column while free water drainage was allowed at the bottom. Once a steady-state velocity profile was established in the sand column, the water supply was stopped (while the air passage was permitted) and the column began to drain from the top due to gravity.
The hydraulic properties of the porous medium were determined by Liakopoulos (1964) for their original experiment and were fitted to and For the gas phase relative permeability, a Brooks-Corey model was used with a pore size distribution index of = 3 , following Geomech. Geophys. Geo-energ. Geo-resour. (2022) 8:107 is constant, gas density is governed by the ideal gas law, where the molar mass of the gas corresponds to air M C = 0.028949 kg mol −1 . Mechanical parameters were not measured by Liakopoulos, therefore we follow the properties suggested by Lewis and Schrefler (1998). Those, and other essential material parameters of the test are given in Table 3.
The numerical model concept is shown in Fig. 8. At top and bottom boundaries, a constant atmospheric gas pressure boundary condition p GR = 101.325 kPa is applied, while capillary pressure at the bottom was fixed to p cap = 0 Pa . Left and right boundaries of the domain are considered impermeable. Displacement in horizontal directions, as well as vertical displacement at the bottom of the column were restricted and a traction free boundary condition is applied at the top of the domain. Initial conditions correspond to fully saturated conditions with atmospheric gas pressure and zero capillary pressure, and the initial stress filed compensates initial pore pressures.
For this test case we apply two solution schemes: a two-phase flow scheme and one with an isobaric gas field (Richards approximation) which has a constant uniform gas pressure of p GR = 101.325 kPa throughout the model domain. The overall simulation time is t = 7200 s , where time step size increases from Δt = 10 s at the beginning to Δt = 400 s at the end of the simulation. The domain is discretized by a one dimensional arrangement of 44 quadrilateral elements.
In Fig. 9, the simulated profiles of gas phase pressure, capillary pressure, saturation, and vertical displacement are shown after t = 5 , 10, 20, 30, 60, and 120 min . At steady-state, the total consolidation settlement at the upper boundary of the column is 1.6 mm while the saturation does not drop below 91 % in both multiphase and Richards variants.
The results of both model variants (two-phase-and Richards flow) agree quite well with those of other authors as the comparison of gas pressure results in Fig. 10 suggests. Another comparison can be found in Fig. 11, where the simulated water pressures of different works are compared with the experimentally determined pressures in the original experiment. Both Figs. 10 and 11 show results from two-phase simulations as well as those using the Richards approximation. As both figures show, the model results of different simulators differ slightly, especially in the transient stage. The reason for this can be found in slightly different model assumptions on which the individual models are built. The differences can be caused by: different minimum permeabilities of the gas phase, use of the ideal gas law or a constant gas density, consideration or omission of mechanical processes, as well as general differences in the numerical formulations (choice of primary variables, coupling schemes, equation solvers, etc.).
It should be noted that it is not the purpose of this article to make a comprehensive comparison of models and to investigate all differences thoroughly. Rather, the comparison of the different solutions is only intended to show that the results of all the works considered (both for the two-phase formulation and for the Richards approximation) show close similarities that confirm the validity of the method proposed here. In this work, we presented a fully coupled, non-isothermal two-phase flow (TH 2 M) model for deformable porous media. Phase transitions via evaporation, condensation and solution processes are considered and follow an equilibrium approach. Governing equations based on fundamental balance equations are developed to weak forms by a standard Bubnov-Galerkin scheme.
The new TH 2 M model is equipped with a comprehensive set of test cases at different coupling levels. Both single processes (T, H, M), binary (TH, TM, H 2 , HM) and ternary (THM, TH 2 , H 2 M) couplings are considered. No test example has yet been created for the complete TH 2 M coupling; however, this will be made up for, either in the form of a (semi)-analytical Four test cases representing a broad range of physical phenomena were selected in this article to verify the implementation of the TH 2 M formulation. These tests were selected because analytical/semi-analytical solutions are available or because they have already been calculated with other codes and thus a basis for code comparison is available. The first benchmark in Sect. 8.1 compared the model results with simple exact solutions for isentropic gas compression to show the validity of the thermodynamic relations. The second test (Point-heatsource, Sect. 8.2) allowed the verification of the non-isothermal hydromechanical process part of the model formulation by comparing the simulation results to an exact solution. The third benchmark test (Heat pipe, Sect. 8.3) examines thermodynamic effects associated with evaporation and condensation of pore water. Finally, the fourth example (Liakopoulos, Sect. 8.4) showed the validity of the two-phase flow formulation and its interaction with the mechanical process. Furthermore, it has been shown that the differences between two-phase flow and the Richards approximation can be significant and have to be taken into account for similar cases, cf. Wang et al. (2011). In addition to the systematic benchmarking procedure, the TH 2 M model has been implemented in an open-source framework (OpenGeoSys) which benefits from professional software engineering (i.e. platform independence, continuous integration, consequent code review for maintaining software quality, automated testing etc.; see Bilke et al. (2019) for details). By choosing an open-source concept for our model, we want to achieve that scientific results can be better understood, acknowledged and improved. The model is freely available along with its test cases, documented and can be modified in order to fit specific needs. In this way, controversial projects, which have a significant impact on the environment, should meet with greater public acceptance.
Further developments of the model presented here are to be carried out in particular for application in the field of nuclear waste disposal. For the investigation of gas transport processes in artificial and natural  phase transitions up to the complete disappearance of one of the phases to be considered. Another aspect of future development is the consideration of enhanced transport in the dilatancy zone induced by very high gas pressures. For this purpose, the mechanical process part is to be extended by damage models and coupled with improved permeability models (c.f. Zill et al. 2021;Yoshioka et al. 2022 Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.

Mass balance without phase transition
Gas pressure equation, where no C-constituent resides in the liquid phase ( C LR = 0 ), thus C GR = GR and x C m,G = 1 , x C m,G = 0.
Capillary pressure equation, where no W-constituent resides in the gas phase ( W GR = 0 ), thus W LR ≡ LR , and molar fractions are always x W m,L = 1 , x C m,L = 0 , x W m,L = 0. (72)