Computational Mechanics (2025) 75:1191–1212 https://doi.org/10.1007/s00466-024-02557-2 ORIG INAL PAPER Phase transition in porous materials: effects of material parameters and deformation regime onmass conservativity Maximilian Brodbeck1 ·Marlon Suditsch1 · Seyed Morteza Seyedpour1,2 · Tim Ricken1,2 Received: 24 March 2024 / Accepted: 13 August 2024 / Published online: 4 December 2024 © The Author(s) 2024 Abstract Phase transition in porous materials is relevant within different engineering applications, such as freezing in saturated soil or pancake sea ice. Mathematical descriptions of such processes can be derived based on Biot’s consolidation theory or the Theory of Porous Media. Depending on parameters such as density ratio, permeability or compressibility of the solid matrix, either small or finite deformations occur. Numerical solution procedures for the general, finite deformation case, suffers from instabilities and high computational costs. Simplifications, assuming small deformations, increases stability and computational efficiency. Within this work shortcomings of simplified theories based on Biot and linearisations of the Theory of Porous Media (TPM) are systematically studied. In order to determine the interaction of the different model parameters a non-dimensional model for poro-elasticity is presented. Based on a characteristic test-case including phase-transition and consolidation, the simplified models are compared to the fully non-linear TPM, focusing on mass errors as well as the time behaviour of the solution. Taking further into account the efficiency of discretisation based on different primal variables and finite-element-spaces, a guideline for selecting an appropriate combination of model, kinematic assumption and discretisation scheme is presented. Keywords Phase transition · Porous material · Homogenisation · Nondimensionalisation · Linearisation · Mass conservation 1 Introduction Fluid flow and deformation are highly coupled in deformable porous materials such as soils, rocks, and tissues. Predict- ing different poro-mechanical processes within small and finite deformation regimes enjoys significant importance in various science and engineering fields, including civil engi- B Tim Ricken tim.ricken@isd.uni-stuttgart.de Maximilian Brodbeck maximilian.brodbeck@isd.uni-stuttgart.de Marlon Suditsch marlon.sudisch@isd.uni-stuttgart.de Seyed Morteza Seyedpour seyed-morteza.seyedpour@isd.uni-stuttgart.de 1 Institute of Structural Mechanics and Dynamics, Faculty of Aerospace Engineering and Geodesy, University of Stuttgart, Pfaffenwaldring 27, 70569 Stuttgart, Germany 2 Porous Media Laboratory, Faculty of Aerospace Engineering and Geodesy, Institute of Structural Mechanics and Dynamics in Aerospace Engineering, University of Stuttgart, Pfaffenwaldring 27, 70569 Stuttgart, Germany neering [1–3], energy and environmental technologies [4–6], and biomechanics [7–10]. Mathematical models for porous materials, which include the conservation laws for mass, momentum, and energy are generally derived through the application of either the continuum theory of mixture (TM) [11–13] respectively the Theory of PorousMedia (TPM) [14, 15], local volume averaging theory [16–18], or Biot’s con- solidation theory [19–22]. The primary distinctions among the models reside in the underlying motivation of the theory and the incorporation of homogenised quantities. Exten- sions of the above mentioned description of porous materials arise, when either phase-exchange or growth affect the phase distribution. Examples for phase-exchange are liquid-to-gas transition in geomechanics [23, 24] or icing (liquid-to-solid transition) in civil and environmental engineering [4, 25–27]. Growth for example occurs within living tissue [7, 28–30]. Within modelling of porous-materials an appropriate dis- cretisation as well as suitable assumptions with respect to the deformation regime are highly relevant. While the Biot the- ory is derived for small deformations, TM aswell as the TPM are applicable also in non-linear regimes. These regimes include geometrically non-linearities (e.g. settlement [31], 123 http://crossmark.crossref.org/dialog/?doi=10.1007/s00466-024-02557-2&domain=pdf 1192 Computational Mechanics (2025) 75:1191–1212 expansion [32], and collapse [33], in geotechnical engineer- ing), plasticmaterial behaviour during shear-bands formation in soils [34] or the descriptions of biological tissue by hyper- elastic material models [9, 35]. Suitable approximations of the deformation state become even more relevant, when cou- pling effects, such as a deformation-dependent permeability [36, 37] or materials with compression point [38] are con- sidered. Beside the kinematic assumptions, discretising the underlying equation system is challenging, due to the saddle point structure of the displacement-pressure system. Stokes stable finite-element spaces are required [39, 40] to avoid spurious pressure oscillations [41, 42]. Classical formula- tions are mainly based on approximations of displacement and pressure, utilising Taylor-Hood elements [43]. Exten- sions to intrinsically unstablefinite-element pairs introducing additional stabilisation can be also found [44]. The accu- racy of such discretisations is impaired due to the violated local mass conservation (relevant e.g. within heterogeneous materials [45–47]), for the description of convection-stable transport processes – which typically required locally con- servative fluxes – or nearly incompressiblematerials [48, 49]. Form an application perspective such robust discretisations are relevant when a large numbers of simulations with model parameters varying over wide ranges are required, such as in polymorphic uncertainty quantification [6, 50] or computa- tional homogenisation of two-scale materials [51]. To retain accurate and computationally efficient solu- tion procedures, suitable combinations of model, kinematic description (small or finite strain theory), and discretisa- tion are highly relevant. The linear Biot theory or at least linerisations of the non-linear TM resp. TPM [26, 42] bene- fits from computationally cheaper and more robust solution algorithms by paying the price due to a limited applicability for finite deformations. This trade-off is the motivation for this paper. Focusing on porous materials undergoing phase- exchange and consolidation, effects of the model parameters on the solution as well as mass conservation are studied. Results gained from the Biot theory, the linerised TPM as well as its fully non-linear formulation are compared. In order to focus on these aspects, the mathematical model is kept simple. Prescribing the phase transition rate, a solution of the energy equation can be avoided. The reduced equa- tion system, based on the TPM, is reviewed in Sect. 2.1. A non-dimensional form, simplified models based on reduced kinematics as well as selected discretisation schemes are presented in Sects. 2.2, 2.3 and 2.4. Utilising these fun- damentals, the intrinsic mass-loss of the simplified models is outlined in Sect. 3. This analysis in then extended to a more realistic case, where phase-transition and consolida- tion occurs simultaneously. Based on a modified Terzaghi problem – the analytic solution is presented in Sect. 4.1 – influences of the non-dimensional model-parameters, a com- parison of simplified and fully non-linear formulations as well as efficiency aspects of selected discretisation methods are presented in Sects. 4.2, 4.3 and 4.4. 2 Material methods 2.1 Governing equations According to the TPM [14, 52] a porous material can be described by a set of superimposed and interacting continua. Fundamentals for the mathematical description, namely the definition of the mixture, its kinematic relations, the relevant balance equations, and the required constitutive models are shortly reviewed within this section. A mixture of immiscible constituents A multi-phasic continuum ϕ is composed of super imposed immiscible continua ϕα . Considering a solid matrix ϕS fully saturated by fluid ϕF, the overall aggregate ϕ is formed by ϕ = ⋃ α ϕα with α = {S, F}. (1) Trough an averaging process, all phase are smeared out over the domain (compare Fig. 1). The ratio between the partial volume dvα of a constituent ϕα and the entire volume dv is described by the volume fraction nα = d vα d v . (2) Based on the mass of a constituent mα , real density ραR and partial density ρα are defined by ραR = dmα d vα and ρα = dmα d v . (3) Their relation can be expressed with respect to the volume fraction ρα = nα ραR. (4) Kinematic relations Each material point P out of a body B is associated with its spatial position x. A set containing all material points is named current configuration �t , while its undeformed state is the reference configuration �0. Field quantities (•) are associatedwith thematerial points. Changes aremeasured by the differential operators gradient and divergence. Gradients are denoted by∇x (•)with respect to current and∇X (•)with respect to to reference configuration, while the divergence is denoted by ∇x · (•) and ∇X · (•). Due to the homogenisation, each position x is associated with two material points PS and PF, with individual motion functions 123 Computational Mechanics (2025) 75:1191–1212 1193 Fig. 1 Averaging process of bi-phasic material. The phases in the real geometry (a) of a representative elementary volume are homogenised (b) and the sum of all partial volumes must be one (c) Fig. 2 Kinematics of a bi-phasic material �α(Pα, t) : { �0 → �t : X(Pα) �→ x (Pα, t ) . (5) This kinematic relation is shown in Fig. 2. Derivation in time leads to the velocity of each constituent ′ xα = ∂�α(Pα, t) ∂t , (6) while the motion of the entire mixture is described by the barycentric velocity ẋ = 1 ρ ∑ α ρα ′ xα. (7) The relative velocity between mixture and individual con- stituent dα = ′ xα − ẋ (8) is named diffusion velocity. Balance equations Balance equations of a mixture, are based on the metaphys- ical principles by Truesdell [53]: 1. All properties of the mixture must be mathematical con- sequences of properties of the constituents. 2. So as to describe the motion of a constituent, we may in imagination isolate it from the rest of the mixture, provided we allow properly for the actions of the other constituent upon it. 3. The motion of the mixture is governed by the same equa- tions as is a single body. The balance equation forψα , a preserved quantity of a phase α, follows (ψα)′α + ψα ∇x · ′ xα = ∇x · φα + σα + ψ̂α, (9) where (•)′α denotes the material time derivative with respect to phase α, ψα the preserved quantity, φα the associated flux, σα the source term and ψ̂ α production due to phase- interaction. Following Truesdells principles, the balance of the mixture is obtained by summation of the phase balances (9). For a closed system this results in ψ̇ + ψ ∇x · ẋ = ∇x · φ + σ . (10) Choosing ψα = ρα results in the balance of mass. For closed systems, the mass of the mixture is preserved, while the mass of each constituent can vary due to phase transition ρ̂α: (ρα)′α + ρα ∇x · ′ xα = ρ̂α. (11) Summation of the component mass balances yields ρ̇ + ρ ∇x · ẋ = 0 , with ∑ α ρ̂α = 0. (12) Assuming incompressible constituents {ρSR, ρFR} = const., (13) the mass balance reduces to (nα)′α + nα ∇x · ′ xα = ρ̂α ραR , (14) which directly determines the volume faction. For ψα = ρα ′ xα , the associated flux can be identified as the Cauchy stress φα = Tα . Considering body forces σα = ραbα the balance of linear momentum yields ρα ′′ xα = ∇x · Tα + ρα bα + p̂α − ρ̂α ′ xα. (15) In analogy to the production of mass, the total production of linear momentum has to vanish: ∑ α p̂α = 0. (16) Assuming incompressible phases (Eq. (13)), creeping pro- cesses ẍ = 0, ′′ xα = 0, ′ xα · ′ xα = 0, (17) 123 1194 Computational Mechanics (2025) 75:1191–1212 and identical body forces acting on each phase bα = b, (18) a equation system for a porous material ∇x · [ TS + TF ] + ρ̂S wFS = s1, ∇x · nF wFS + ∇x · ′ xS = s2, (nS)′S + nS ∇x · ′ xS = s3. (19) Herein, (19)1 and (19)2 are the mass and momentum balance of the entire mixture and (19)3 the mass balance of the solid phase. The source terms are defined by s1 = −ρ b, s2 = s3 [ 1 − ρSR ρFR ] , s3 = ρ̂S ρSR , (20) where ρ denotes the density of the mixture ρ = ρFR + nS [ ρSR − ρFR ] . (21) Constitutive relations In addition to the balance Eq. (19), constitutive relations for the Cauchy stresses TS and TF, the volume flux nF wFS and the phase transition rate ρ̂S are required. In order to ensure thermodynamically consistent solutions, these are evaluated based on the Clausius-Duhem inequality. Under the assump- tions Eqs. (13), (17) and (18), the inequality simplifies to ∑ α [Tα E · Lα − ρα (ψα)′α − ρ̂α ψα] − p̂FE · wFS ≥ 0, (22) with TS E = TS + nS p I, TF E = TF + nF p I, p̂FE = p̂F − ρ̂α ′ xα − p∇x nF. , (23) where p denotes the fluid pressure. According to Simo and Pister [54] and Diebels et al. [55] the extra stresses follow from the reversible part of the Clausius-Duhem inequality (22) by TS E = μS [ I − C−1 S ] + λS ln(JS) and TF E = 0, (24) wherein μS and λS denote the first and second Lamé param- eter of the solid phase. The remaining dissipative parts of the Clausius-Duhem inequality (22) D = Dρ̂ + Dp̂ ≥ 0, (25) are used to determine volume flux and transition rate. For constant transition rates,Dρ̂ > 0 is fulfilled [56]. Postulating a linear relation between p̂FE and seepage velocity p̂FE = −μFR kS0S wFS, (26) Dp̂ ≥ 0 and thus the dissipation inequality (25) is valid. The isotropic model parameters are the intrinsic permeability kS0S and the dynamic viscosity of the fluid μFR. Inserting p̂FE into the momentum balance of the fluid (15) yields the volume flux nFwFS = − kS0S μFR ( ∇x p + ρFR b − ρ̂S nF ′ xS ) . (27) Within the following, kS0S/μ FR, the ratio of intrinsic perme- ability andfluid viscosity, is theDarcy permeability kD,while the body force b purely depends on the gravity constant g. 2.2 Nondimensionalisation Following Buckinghams -Theorem a set of equations, con- taining n parameters with rank k of the resulting dimension matrix, can be expressed with respect to n = n − k (28) dimensionless parameters i . Non-dimensional parameters for poro-elastic problems can be found in Terzaghi [57], Sondberg [58] or Osman and Randolph [59]. Extending their work to problems containing mass exchange – the model parameters are summarised in Table 1 – results in n = 10 non-dimensional parameters: x̃ = x � , , t̃ = t kD μS �2 , ũS = uS � , p̃ = p μS , 1 = λS μS , 2 = ρSR ρFR , 3 = ρFR � g μS , b̃ = b g , 4 = �2 μS ρFR ( kD )2 , 5 = kD ρ̂S . (29) In order to derive the non-dimensional form of the TPM, the dimensionless parameters (29) are calculated based on reference values (•) ∣∣ re f for each material parameter. 123 Computational Mechanics (2025) 75:1191–1212 1195 Table 1 Dimension table for a poro-elastic model containing mass exchange x t uS p ρSR λS μS ρFR kD = kS0S μFR b ρ̂S g � M 0 0 0 1 1 1 1 1 -1 0 1 0 0 L 1 0 1 -1 -3 -1 -1 -3 3 1 -3 1 1 T 0 1 0 -2 0 -2 -2 0 1 -2 -1 -2 0 Defining their non-dimensional form ˜(•) = (•) (•) ∣∣ re f , (30) and transforming spatial- andmaterial time derivatives based on ∂(•) ∂x = 1 � ∂(•) ∂ x̃ , (•)′S = kD ∣∣ re f μS ∣∣ re f �2 〈•〉′S , (31) the non-dimensional equation system of the TPM reads ∇x̃ · [ T̃S + T̃F ] + 5 w̃FS = s̃1, ∇x̃ · 〈ũ〉′S + ∇x̃ · nFw̃FS = s̃2, 〈 nS 〉′ S + nS ∇x̃ · 〈ũ〉′S = s̃3. (32) The non-dimensional source terms are defined by s̃1 = − 3 ρ̃ b̃, s̃2 = [ 1 ρ̃SR − 2 1 ρ̃FR ] 4 5 2 , s̃3 = 1 ρ̃SR 4 5 2 , (33) where ρ̃ denotes the non-dimensional density of the mixture ρ̃ = ρ̃FR + nS ( 2 ρ̃SR − ρ̃FR ) . (34) Analogously, the constitutive relations can be transformed. The Cauchy stresses (23)1,2 yields T̃S + T̃F = μ̃S [ I − C̃−1 S ] + 1λ̃ S ln ( J̃S ) I − p̃ I (35) while the seepage velocity (27) results in nFw̃FS = −k̃D ∇x̃ p̃ + nFw̃ext FS , nFw̃ext FS = k̃D [ 3 ρ̃FR b̃ − 5 nF 〈ũS〉′S ] . (36) For constant material parameters the governing equations can be further simplified. Choosing (•) = (•) ∣∣ re f , the non- dimensional material parameters λ̃S, μ̃S, k̃D, ρ̃SR, and ρ̃FR are equal to one. 2.3 Pull-back and linearisation In the previous sections, the governing equations of a porous material based on the TPM were derived and non- dimensionalised. For a solution using the finite element method (FEM) and a total Lagrangian ansatz, the gov- erning equations in current configuration (32) have to be pulled back into reference configuration. Solving this system of non-linear equations is computationally complex. Under the assumption of small deformations, simplifications are possible. In this section, descriptions based on the fully non- linear TPM (37), linerarisations followingRicken (lTPM(R)) respectively Bluhm (lTPM(B)) (45), and Biot (42) are pre- sented. The fully non-liner case Applying the required transformations on the governing equations in current configuration (32) leads to ∇X̃ · [ P̃S + P̃F ] + J̃S 5 w̃FS = S1, ∇x̃ · 〈ũ〉′S + ∇X̃ · nF w̃FS,0S = S2, 〈 nS 〉′ S + nS J̃S ∇x̃ · 〈ũ〉′S = S3. (37) The source terms in reference configuration can be written as S1 = J̃S s̃1, S2 = J̃S s̃2, S3 = s̃3. (38) Mapping stresses, the divergence of the solid velocity, and the seepage-velocity from current- into reference configuration yields P̃S + P̃F = J̃S [ T̃S + T̃F ] F̃ −T S , (39) ∇x̃ · 〈ũ〉′S = J̃S ∇X̃ 〈ũ〉′S · F̃−T S , (40) nFw̃FS,0S = J̃S [ −k̃D GP + nFw̃ext FS ] , (41) 123 1196 Computational Mechanics (2025) 75:1191–1212 where P̃S + P̃F denotes the first Piola-Kirchhoff stress of the mixture. For better readability the abbreviation GP = ∇X̃ p̃ · C̃−1 S is introduced. Linearisation for small deformations In situations, where only small deformations are present, the equation system (37) can be simplified. Before discussing linearisations of the TPM, an extension of the Biot equations is presented as a simple alternative. The Biot equations can be recasted from the TPM, assum- ing a linear stress–strain relation, F̃S = I, and J̃S = 1: ∇X̃ · ( T̃S + T̃F )B + 5 (w̃FS) B = s̃1, ∇X̃ · 〈ũS〉′S + ∇X̃ · ( nFw̃FS )B = s̃2, 〈 nS 〉′ S + nS ∇X̃ · 〈ũS〉′S = s̃3. (42) The source terms followEq. (33),while the linear constitutive relations for the second Piola-Kirchoff stress of the mixture and the seepage velocity are defined by: ( T̃S + T̃F )B = 2μ̃S ε̃S + 1 λ̃S ∇X̃ · ũS I − p̃ I, (43) ( nFw̃FS )B = −k̃D ∇X̃ p̃ + nFw̃ext FS . (44) Linearisations of the fully non-linear TPM (37) are pro- posed byRicken [42] (lTPM(R)) andBluhm [26] (lTPM(B)). The approach by Ricken is based on a systematic linearisa- tion, calculating a Taylor series, which is truncated after the linear element around the reference configuration. In con- trast, Bluhm proceeds heuristically, replacing each term by it’s linearised counterpart. This preserves the structure of the non-linear equation system, but has no clear theoretical foun- dation. Main difference between lTPM and Biot results from a formal lineraisation of the mapping between the configu- rations. This is entirely neglected within the Biot equations, where the alikeness of current and reference configuration is assumed. For a better readability the following abbreviations are introduced: P = LIN [ P̃S + P̃F ] , Q = LIN [ J̃S 5 w̃FS ] , V = LIN [∇x̃ · 〈ũ〉′S ] , W = LIN [ nF w̃FS,0S ] , J = LIN [ J̃S ] , ε̃S = sym [∇X̃ ũS ] . Thegeneral structure of the linearised equations then follows: ∇X̃ · P + 5 Q = J s̃1, V + ∇X̃ · W = J s̃2, 〈 nS 〉′ S + nS ∇X̃ · 〈ũS〉′S = s̃3. (45) Effects of the different linerisations on J , P , Q, V and W are summarised in Table 2. 2.4 Weak forms and numerical treatment Solving a partial differential equation (PDE) system with the FEM requires the weak form of the governing equa- tion system. Primal variables are approximated in finite- element-spaces (FE-spaces). Suitable Inf-Sub conditions avoid numerical instabilities. Throughout this analysis dis- cretisations based on • displacement, (fluid) pressure and volume fraction (upn- formulation) • displacement, (fluid) pressure, total pressure and volume fraction (uppn-formulation) • displacement, seepage velocity, (fluid) pressure, and vol- ume fraction (uvpn-formulation) are considered. Within this section the resulting weak forms, as well as requirements on the FE-spaces are discussed. Time derivatives are thereby generally captured using backward differentiation formulas of order 1 (BDF1). The upn-formulation Choosing displacement, pressure and volume fraction as pri- mal unknowns, requires a solution (ũS, p̃, nS) ∈ ( H1 )d × H1 × L2 such that a1([ũS, p̃], vu) + a2([ũS, nS], vu) = l1(vu), b1(ũS, vp) + b2([ũS, p̃], vp) = l̂2(vp), c1([ũS, ũnS, nS], vn) + c2(nS, vn) = l̂3(vn), (46) holds for all test function (vu, vp, vn) ∈ ( H1 )d × H1 × L2. The modified right-hand-sides l̂ i , required due to the time discretisation, are defined by l̂2(vp) = ∫ �0 ∇x̃ · 〈ũ〉′S · vp dV + l2(vp) , (47) l̂3(vn) = c2 (( nS )n , vn ) , (48) where ũnS and ( nS )n denotes solutions of the last converged time-step. Let now �0 be the reference configuration and � ⊆ ∂�0, the forms ai , bi , ci , and li are defined by a1 = ∫ �0 [ P̃S + P̃F ] · ∇X̃ vu dV , a2 = ∫ �0 − [ J̃S 5 w̃FS − S1 ] · vu dV , l1 = ∫ � [ P̃S + P̃F ] · N · vu dA , (49) 123 Computational Mechanics (2025) 75:1191–1212 1197 Table 2 Linerised TPM models following Ricken (lTPM(R)) and Bluhm (lTPM(B)) lTPM(R) lTPM(B) J 1 + ∇X̃ · ũS P 2μ̃S ε̃S + 1 λ̃S ∇X̃ · ũS I − h1 p̃ Q − k̃D nF [ h2 · ∇X̃ p̃ + J 3 ρ̃FR b̃ − 5 nF 〈ũS〉′S ] W nF V ∇X̃ · 〈ũS〉′S J ∇X̃ · 〈ũS〉′S W −k̃D [ ∇X̃ p̃ · h1 + J 3 ρ̃FR b̃ − 5 nF 〈ũS〉′S ] −k̃D∇X̃ p̃ · h1 + J nFw̃ext FS h1 J I − 2 ε̃S J I h2 ∇X̃ ũS − J I - b1 = ∫ �0 [∇x̃ · 〈ũ〉′S − S2 ] · vp dV , b2 = ∫ �0 − t̃ nFw̃FS,0S · ∇X̃ vp dV , l2 = ∫ � − t̃ nFw̃FS,0S · N · vp dA , (50) c1 = ∫ �0 [ nS F̃ −T S · F̃S − t̃ S3 ] · vn dV , c2 = ∫ �0 nS · vn dV , (51) using F̃S = ∇X̃ [ ũS − ũnS ] . The primal fields are typically approximated using Lagrangian elements. While equal order (EO) approximations for displacement and pressure are not Stokes stable, Taylor-Hood (TH) or Mini elements are more suitable. The uppn-formulation In order to improve robustness of a discretisation in the incompressible limit, an additional total pressure p̃t = p̃ − 1λ̃ S ∇X̃ · ũS (52) can be introduced [48, 49]. The weak forms of this method, which is only applicable for theBiot equations (42), are given in Appendix 1. Beside the parameter robust a-priori estimate, this scheme decouples the FE-spaces of displacement and fluid pressure. Stokes stability is only required for the spaces of displacement and total pressure, e. g. TH ormini elements. This enables stable and equivalently accurate approximations for displacement and pressure of higher order. The uvpn-formulation The so far discussed methods are not locally mass conserva- tive. This can be achieved by introducing the seepage velocity as additional primal variable. Local mass conservativity can be relevant when discontinuities e. g. in permeability or phase transition occurs. Within application to the (l)TPM a direct discetisation of the term J̃S 5 w̃FS is not pos- sible. As this term is mostly neglected in literature [60, 61], this seems acceptable. The weak forms of this method Fig. 3 A 2D unit square (thickness h, symmetry boundary conditions) with prevented outflow are given in Appendix 1. FE-spaces are typically based on Lagrangian elements for the displacement and lowest order Raviart-Thomas (RT) elements for the flow alongside with element-wise constant pressure approximations [62]. Higher orders are possible, but require non-conforming displace- ment fields [63, 64]. 3 A first glance onmass losses in linearised poro-elasticity To gain first insights into the problem of mass losses a unit square (thickness h) under symmetry boundary conditions with prevented outflow (Fig. 3) is analysed. For constant tran- sition rate ρ̂S, the displacement varies linearly, while the pressure and volume fraction are constant in space. The current position of a material point xS can be expressed as a linear function of its reference configuration XS xS = [ (1 + ã) XS,x (1 + ã) XS,y ] , (53) where ã denotes the equal elongation in the two space dimen- sions, normalised by the domain-length �. Evaluating the overall mass balance of each model (37)2, (45)2 resp. (42)2, 123 1198 Computational Mechanics (2025) 75:1191–1212 Table 3 The reformulated overall mass balances as well as the dis- placement solutions for the 2Dphase-transition-problemusing different models Model total mass balance ã TPM 2ã′ (1 + ã) = (1 + ã)2 c exp ( c t̃/2 ) − 1 lTPM(R) 2ã′ (1 + ã) = (1 + 2ã) c [ exp ( c t̃ ) − 1 ] /2 lTPM(B) 2ã′ = c c t̃/2 Biot 2ã′ = c c t̃/2 results in an ordinary differential equation (ODE) for the elongation a. These ODEs as well as their solutions for the initial-value ũS ( X̃S, t̃ = 0 ) = 0 resp. ã ( t̃ = 0 ) = 0 are summarised in Table 3. For the sake of clarity the factor c = (1 − 2) 4 5 −1 2 (54) is introduced. Inserting the displacement solution (see Table 3) into the solid mass balance of each model (37)3, (45)3 resp. (42)3 yields 〈 nS 〉′ S + nS c − 4 5 2 = 0. (55) This model-independent ODE results under the initial con- dition nS(X̃S, t̃ = 0) = nS0S to nS = { c1 exp(−c t̃) + c2, 2 �= 1 nS0S + 4 5 ( 2) −1 t̃, else , (56) with c1 = nS0S − c2 and c2 = (1 − 2) −1. Amore detailed derivation of these ODEs for the case of the fully non-linear TPM can be found in the Appendix 1. Based on the non-dimensional mass M̃ = M h�2 ρFR = (1 + ã)2 [ 1 + ( 2 − 1) nS ] , (57) a relative mass error is defined by errM = M̃(t̃) − M̃(t̃ = 0) M̃(t̃ = 0) . (58) Figure4 shows the mass error of the different models. It is obvious, that conservation is only fulfilled by the TPM. The mass error of the small deformation models results purely from the insufficient approximation of the volume dilation, as the volume fraction is independent of the model. Fig. 4 Mass error errM within in a unit square undergoing phase tran- sition, modelled by different sets of governing equations 4 Results and discussion In the previous section a mass error due to insufficient approximations of the volume dilatation during phase tran- sition was introduced. The analysis is now extended to more general cases, by incorporating fluid flow. In this context the following hypothesis is formulated: The interaction of phase transition and consolidation is typically dominated by one of the two processes. Dif- ferent parameter combinations favour the domination of either phase transition or consolidation and affect the suitability of small deformation models. In order to examine this hypothesis, the following questions have to be investigated: 1. Howmodel parameters affect displacements and therefore the assumption of small deformation? 2. How does the different models for small deformation affect the dynamic solution behaviour compared to the fully non-linear TPM? To answer these questions, a variant of Terzaghi’s consol- idation problem [65] is studied, from which the traction boundary condition is removed while phase transition is incorporated. The derivation of an analytic solution, pre- sented in Sect. 4.1, allows the systematic determination of influences of the model parameters, presented in Sect. 4.2. The results of the different models at different deformations stages are compared in Sect. 4.3. As practice relevant prob- lems are typically computationally expensive, accuracy and efficiency of establishes, as well as currently developed dis- cretisations schemes are compared in Sect. 4.4 and closes the 123 Computational Mechanics (2025) 75:1191–1212 1199 Fig. 5 Geometry and boundary conditions of amodified Terzaghi prob- lem. The traction condition is removed, while phase transition is applied examination of the hypothesis. All numerical results within this section were produced using DOLFINx [66], while the arising linear equation systems are solved using the direct solver MUMPS [67]. 4.1 Problem and analytic solution The basis for the following analysis is a slender domain (�/L = 10) with spacial- and temporal constant phase transition. The boundary conditions are depicted in Fig. 5. Due to the applied conditions, the solution simplifies into a quasi one-dimensional (1D). Additionally it is assumed, that 5 w̃FS resp. ( 5/n F) 〈ũS〉′S can be neglected. Clearly such a solution can only be found for the linearBiot equations (42). The governing equations within 1D reduce to ( 1 + 2) ∂2ũ ∂ (z̃)2 − ∂ p̃ ∂ z̃ = 0, ∂2ũ ∂ t̃ ∂ z̃ − ∂2p̃ ∂ (z̃)2 − (1 − 2) 4 5(t) 2 = 0, ∂nS ∂ t̃ + nS ∂2ũ ∂ t̃ ∂ z̃ − 4 5(t) 2 = 0. (59) Based on the idea of Stickle and Pastor [65], solutions for displacement and pressure can be obtained by integrat- ing the balance of linear momentum (59)1 over z̃, taking its time derivative and inserting the result into the overall mass balance of the mixture (59)2. The resulting Poisson- like equation ∂ p̃ ∂ t̃ = ( 1 + 2) ∂2p̃ ∂ (z̃)2 + (1 − 2) 4 5(t) 2 (60) can be solved under the given initial- and boundary condi- tions in two steps: For the homogeneous case, the solution of the underlying Sturm-Liouville problem follows p̃(z̃, t̃) = ∞∑ n=1 pn(t̃) yn(z̃), (61) where yn are time-constant eigenfunctions yn = sin ( (1 + 2n)π (1 − z̃) 2 ) , (62) while the eigenvalues pn are time-dependent. Inserting Eqs. (61) into (60) and expressing the phase transition term as Fourier series based on the eigenfunctions (62) leads to an ODE ac:ODE in time for the n-th eigenvalues of the pressure ∂ pn ∂ t̃ + Nn pn = (1 − 2) 4 5(t) 2 , (63) where Nn is defined by Nn = ( 1 + 2) [ (1 + 2n)π 2 ] . (64) Assuming a spatially constant phase transition term varying in time 5(t) = { 5 t̃ ≤ t̃Gr 0 else , (65) Eq. (63) can be solved by pn(t) = 2 [ 2 (1 + 2n)π ]3 c ψn(t), (66) introducing a time function ψn defined by ψn = { 1 − exp (−Nnt) t ≤ tGr exp ( τ̃Gr,n ) − exp ( t̃n ) else t̃n = −Nn t̃, τ̃Gr,n = t̃n + Nn t̃PT. (67) The solution of the displacement can be calculated from Eq. (59)1 by integrating the pressure solution in space: ũ(z̃, t̃) = 1 1 + 2 ∞∑ n=1 2pn(t̃) (1 + 2n)π ŷn(z̃), ŷn = cos ( (1 + 2n)π (1 − z̃) 2 ) . (68) As a last step the solid volume fraction can be calculated based on displacement and pressure. Up to the authors best knowledge there is no applicable analytic solution, due to 123 1200 Computational Mechanics (2025) 75:1191–1212 the combination of non-homogeneity and non-constant coef- ficient within the solid mass balance (59)3. Therefore an numeric ODE solver is applied, to calculate the volume frac- tion. Determination of mass and mass-loss In order to evaluate mass, outflow and the overall mass error, the determinant of the deformation gradient J̃S as well as the inverse of the right Cauchy-Green deformation tensor C−1 S are required. In the considered 1D case, these are solely dependent on the the derivative of the ũ with respect to z̃: ∂ ũ ∂ z̃ = p̃ 1 + 2 . (69) Based thereon, as well as on the solid volume fraction nS, the mass within the domain M̃ = M A� ρFR = ∫ 1 0 J̃S [ 1 + ( 2 − 1) nS ] dz̃, (70) and the outflux, pulled back into reference configuration (41), nFw̃FS,0S = − [ 1 + ∂ ũS ∂ z̃ ]−1 ∂ p̃ ∂ z̃ (71) are evaluated. Integrating the outflux over surface and time yields the outflown mass M̃out(t̃) = ∫ t̃ 0 nFw̃FS,0S(�, τ ) · N dτ = ∫ t̃ 0 −∂ p̃ ∂ z̃ (�, τ ) dτ, (72) where A is the constant cross-section of the 1D problem and N its normal vector in reference configuration. Having these definitions at hand, the relative mass error (58) is defined by errM = M̃(t̃) + M̃out(t̃) − M̃(t̃ = 0) M̃(t̃ = 0) . (73) 4.2 Influence of model parameters on themass loss Based on the analytic solution, the significance of the mass error depending on the combination of the non-dimensional parameters has to be evaluated. Within eqs. (59), the com- pressebility paramter 1, the density ratio 2 and the phase-transition-rate 4 5 can be varied independently. For a better clarification of the effects, detached from physically motivated parameter spaces, the ranges, specified in Table 4, are considered. The density ratio 2 is limited to values smaller than 1 to avoid situations, where fluid is soaked into the domain and then converted into solid. Table 4 Parameter ranges of the varied non-dimensional parameters 1, 2 and 4 5 1 2 4 5 min max min max min max 0.25 10 0.2 0.9 0.01 10 Fig. 6 Maximum displacement at t̃ = t̃PT for 2 = 0.2. Phase transi- tion ( 4 5) and compressibility ( 1) are varied To make the different cases comparable, the fully con- solidated state is evaluated while the transition time t̃PT is adjusted, such that the solid mass is identical. Therefore, an analytic estimate based on the initial mass M̃(t̃ = 0) and the outflow (72) is used: M̃(t̃ → ∞) = M̃(t̃ = 0) − M̃out(t̃ → ∞) ↔ M̃S(t̃ → ∞) = nS0S � + 4 5 t̃PT 2 . (74) Prescribing the ratio between final- and initial solid mass, the transition time can be estimated by t̃PT = nS0S � 2 4 5 [ M̃(t̃ → ∞) M̃(t̃ = 0) − 1 ] . (75) Within all following calculations nS0S = 0.2 and M̃S(t̃ → ∞) = 0.75 A� is assumed. Displacement and mass error for the lowest density ratio ( 2 = 0.2) of this study are displayed in Figs. 6 and 7. In order to additionally determine the effect of the density ratio 2, two bounding curves of the discussed parameter space are selected. The first bounding curve is defined by a fixed compressibility parameter ( 1 = 0.25) and a varying transition rate 4 5. In contrast, the second one is based on 4 5 = 10, while 1 is varied. Effects on the mass error are shown in Fig. 8 respectively 9. It becomes obvious that high transition rates ( 4 5) for compressible materials (low 123 Computational Mechanics (2025) 75:1191–1212 1201 Fig. 7 Mass error in the fully consolidated state for 2 = 0.2. Phase transition ( 4 5) and compressibility ( 1) are varied Fig. 8 Mass error in the fully consolidated state for 1 = 0.25. Phase transition ( 4 5) and density ratio ( 2) are varied 1) with a significant density ration (low 2) leads to high mass errors. This coincides with large displacements directly after the phase transition ends (exemplarily shown in Fig. 6). The presented results show clearly, that the linear theory after Biot has an inherentmass error, when outflow is present. In order to determine the influence of the flux approxima- tion, the outflow boundary condition is removed from the top surface of the domain, and displacement and pressure are evaluated analogously as in Sect. 3. The relative mass error (58) then yields errM, NO = exp (−c t̃PT ) [ c t̃PT + 1 ] − 1. (76) A direct comparison between the phase transition prob- lem with and without outflow, keeping the non-dimensional parameters as well as the phase transition time identical, is shown within Figs. 10 - 12. Therein the difference of the mass error errM = errM −errM, NO is depicted. Phase tran- Fig. 9 Mass error in the fully consolidated state for 4 5 = 10. Com- pressibility ( 1) and density ratio ( 2) are varied sition and compressibility are varied for a fixed density ratio 2 within Fig. 10. The two bounding curves, varying 2, are shown in Figs. 11 and 12. The results show three differ- ent regions. With increasing transition-rate, phase transition dominates over consolidationwhich results in large outfluxes at relatively high displacements. Even if these displacements are significantly smaller in cases with outflow, the overall error is larger. This behaviour indicates that, in addition to the wrong approximation of the volume dilatation, there is an additional error due to the wrong flux approximation. This error is relevant, whenever large deformations and signifi- cant outflows occur. A second region can be identified for small transition ratios, which results in large dilation errors for cases without outflow. The error for the outflow cases becomes smaller than for the case without outflow. This results from the rapid consolidation, preventing large defor- mations and therefore larger errors. Thirdly and finally nearly no or if so a slightly larger error in cases without outflow can be seen formoderate density ratios near one.Major limitation within this comparison are the quite different displacement state between the cases with identical parameters. Within a second comparison the transition time t̃PT is adjusted for the case without outflow, such that the displacement at t̃ = t̃PT is identical in the case with and without outflow. It becomes apparent that the flow-dependent error dominates the cases with outflow boundary condition. 4.3 Influence of themodel on the solution Up to this point only the Biot equations were analysed.While most of the results regarding accuracy and efficiency can be transferred from the Biot equations to the (l)TPM – only the four field formulations based on displacement, pressure, total pressure and volume fraction can not be transferred directly – main focus of this section is an analysis of influences of the 123 1202 Computational Mechanics (2025) 75:1191–1212 Fig. 10 Difference in mass error for problem with and without outflow for 2 = 0.2. Phase transition ( 4 5) and compressibility ( 1) are varied Fig. 11 Difference in mass error for problem with and without outflow for 1 = 0.25. Phase transition ( 4 5) and density ratio ( 2) are varied different models. The considered test-cases, based on 1 = 0.9 and 2 = 0.7, are summarised in Table 5. For the cases 1-x the transition rate 4 5 is varied, estimating t̃PT by Eq. (75), while for the cases 2-x the transition rate 4 5 is fixed to 10. Thereby, the transition time is adjusted, such that a specified maximal displacement occurs. All cases are solved on a mesh of 6 × 160 cells, each containing two triangular elements. A classical upn-formulation and a Taylor-Hood pair of order one for displacement and pressure is used. The volume-fraction is approximated using Lagrangian elements of order one. The constant time-step t̃ is chosen, such that Fig. 12 Difference in mass error for problem with and without outflow for 4 5 = 10. Compressibility ( 1) and density ratio ( 2) are varied Table 5 Testcases for the model-comparison with maximal displace- ments ũmax between 0.01 and 0.1 id 4 5 t̃PT t̃ ũmax 1-1 0.21 1.83 5.09 · 10−5 0.01 1-2 0.63 0.61 2.75 · 10−5 0.03 1-3 1.2 0.32 5.65 · 10−5 0.05 1-4 4.1 0.094 7.57 · 10−5 0.1 2-1 10 0.0025 8.33 · 10−5 0.01 2-2 10 0.0079 8.78 · 10−5 0.03 2-3 10 0.014 9.20 · 10−5 0.05 2-4 10 0.03 9.38 · 10−5 0.1 the solution using the Biot equations fulfils ∣∣∣ũS,ext − ũhS ∣∣∣ 1,0 < 10−3, ∣∣∣p̃ext − p̃h ∣∣∣ 1,0 < 1.5 · 10−2, ∣∣∣∣ ( nS ) ext − ( nS )h∣∣∣∣ 0 < 10−3. To determine the influence of the models on the solution, the maximal displacement (nondimensionalised) ũmax = ũ(�, t̃) as well pressure p̃max = p̃(0, t̃), the systemmass M̃ as well as the mass error are evaluated. The evolution of ũmax and p̃max are shown in Fig. 13. Comparing the four cases, which have by intention a comparably small mass error, the most notable difference is the evolution of displacement and pressure especially for the cases 1-2 and 1-3. Focusing on the displacement (Fig. 13 left column) and comparing ũmax, up to which the simplified solutions show a similar behaviour than 123 Computational Mechanics (2025) 75:1191–1212 1203 Fig. 13 ũmax (left) and p̃max (right) for the cases 1-n (row n) calculated using Biot and the (l)TPM 123 1204 Computational Mechanics (2025) 75:1191–1212 Fig. 14 ũmax (left) and p̃max (right) for the cases 2-2 (row 1) and 2-3 (row 2) the non-linear TPM, the maximum displacement is signifi- cantly higher for the case 1-4 that for 1-2 resp. 1-3. Similar observations can be made when p̃max is compared, even though the pressures start to diverge earlier, than the dis- placement. This observation can even be supported, when, instead of adjusting the transition parameter 4 5, the tran- sition time is adjusted, while keeping 4 5 = 10 constant (test-cases 2-x). Results for the cases 2-2 and 2-3 are depicted within Fig. 14, highlighting the fact, that for equal ũmax, the simplified solutions match the TPM quite accurately. The fluid pressure is overestimated but at least the shape of the evolution in time is captured quite accurately. These findings are somehow a contradiction to small deformation models in classical elasticity. While they suite non-linear theories typically up to a certain strain level, the behaviour between phase transition and consolidation seems to be important within poro-elasticity with phase transition. While within case 1-1 and 1-4 consolidation is either very fast or very slow compared to phase transition, the cases 1-2 and 1-3 are in some sense a superimposition. The serious underestimation of the consolidation time of the simplified models, compared to theTPM, causes significant differences, not onlywith respect to extremal values, but alsowith respect to the temporal evolution. Up to this point displacement and pressure were analysed, while neglecting influences on the system mass. Due to the results within Sect. 3, mass errors for the linearised TPM just like for the Biot problem are expectable. Mass errors for Biot, as well as the linearised TPM by Bluhm (lTPM(B)) are given in Table 6. For this 1D cases, the linearised TPM following Ricken (lTPM(R)) is free of mass losses. Based on the findings from Sect. 3, this will not hold for 2D resp. 3D cases. The same observation holds true in the limit case ( 1 = 0.25, 2 = 0.2 and 4 5 = 10). While no mass error occurs for calculations based on the fully non-linear TPM as well as lTPM(R), 6.9% resp. 6.2% of the initial 123 Computational Mechanics (2025) 75:1191–1212 1205 Table 6 Negative mass error −errM in [%] at t̃ = t̃PT id Biot lTPM(B) id Biot lTPM(B) 1-1 1.7E-1 1.3E-2 2-1 5.7E-3 5.4E-3 1-2 4.2E-1 6.6E-2 2-2 5.2E-2 4.7E-2 1-3 6.1E-1 1.8E-1 2-3 1.5E-1 1.3E-1 1-4 9.0E-1 5.7E-1 2-4 6.5E-1 5.3E-1 mass are lost, using Biot resp. lTPM(B). Focusing beside this actual mass loss on the evolution of the system mass, the overestimated consolidation becomes one again visible. Figure15 shows the entire mass within the domain over the simulation time. While nearly no difference for the total mass is present in case 1-1, a significant underestimation of the mass, especially during consolidation, is observed in case 1-4. Similar tendencies can be observed for parameter combinations with higher mass errors. The rapidity of the consolidation is thereby increasingly overestimated, when larger displacements occurs. 4.4 Efficiency of numerical solution procedures In the previous sections an analytic solution for the 1D phase transition problem, using the Biot equations, as well as influences of model parameters on the solution resp. the inherentmass errorwere discussed. In amore general setting, either with respect to boundary conditions or when non- linearities arise within the governing equations, numerical solution procedures are required. Therefore the efficiency of the different solution procedures, discussed in Sect. 2.4, will be evaluated. For this convergence study, the bound- ary value problem (BVP) (see Fig. 5) with model parameters 1 = 0.25, 2 = 0.2 and 4 5 = 10 and a phase transition time t̃PT = 0.011 is considered. The domain is discretised using triangular elements in space and a fixed increment t̃ = t̃PT/320 in time. A summary of considered discretisa- tion schemes is given in Table 7, the sequence of considered meshes as well as the number of degrees of freedom nDOF in Table 8. Fig. 15 The total mass M̃ for the cases 1-1 (top left), 1-2 (top right), 1-3 (bottom left) and 1-4 (bottom right) calculated using Biot and the (l)TPM 123 1206 Computational Mechanics (2025) 75:1191–1212 Table 7 Overview over the tested discretisations, summarising abbre- viation (ID), governing equations (Gov. Equation) and considered finite-element spaces (FE spaces) ID Gov. Equation FE spaces upn-mini-0 (46) (P1 + B)2 × P1 × DP0 upn-mini-1 (46) (P1 + B)2 × P1 × (D)P1 upn-TH1 (46) (P2)2 × P1 × (D)P1 uvpn-1 (A9) (P1)2 × RT0 × DP0 × DP0 uvpn-2 (A9) (P2)2 × RT0 × DP0 × (D)P1 upn-EO2 (46) (P2)2 × P2 × (D)P1 upn-TH2 (46) (P3)2 × P2 × (D)P2 uppn-TH1-2 (A1) (P2)2 × P2 × P1 × (D)P1 uppn-TH2-3 (A1) (P3)2 × P3 × P2 × (D)P2 For a systematic comparison of the different discretisation schemes, norms of the displacement, pressure and volume fraction error are evaluated at t̃ = t̃end. Except for the uvpn- 1 and uvpn-2 formulation, displacement and pressure error are measured within the H1-seminorm | f |1,0 = ∫ �0 ∇X f · ∇X f dV , (77) while the error of the volume fraction and the seepage veloc- ity nFw̃FS,0S are evaluated within the L2-norm | f |0 = ∫ �0 f · f dV . (78) For the uvpn-1 and uvpn-2 formulation, the displacement error is evaluated in the H1-seminorm, while the errors of pressure, volume fraction and seepage velocity are evaluated in the L2-norm. The convergence history of the field quantities calcu- lated by the different schemes is depicted in Fig. 16. The errors are plotted over the solution time per newton iteration, which contains the time required for reassemblying the tan- gent and solving the linear equation system. The presented comparison is based on solution times using discontinuous approximations for the volume fraction. Even if the number of degrees of freedom is higher for these schemes, the overall solution time is at least equal or smaller than the comparable, continuous approximation. The actual error is for the here considered case identical for continuous and discontinuous approximation. Focusing at first on the spacial discretisation, quite well known shortcoming can be shown for the upn-based for- mulations. If a Taylor-Hood pair is used for displacement and pressure, the convergence ratio of the pressure is one order lower than of the displacement, requiring fine meshes which are for the approximation of displacement or volume fraction not necessary. This effect is stronger for the low- est order case (upn-TH1) and alleviate for the higher order case (upn-TH2). Improvements are possible, if equal approx- imations orders are used for displacement and pressure. To achieved first order accuracy for displacement and pressure, the Taylor-Hood pair can be be replaced by a Mini element. The comparably cheap displacement approximations pays of, when the pressure is solved down to a specified toler- ance. Comparing the upn-mini-based formulations with the upn-TH1 formulation, the solution time is between 40 an 50% smaller. Formulations using equal, but higher orders (e.g. upn-EO2, uppt-TH1-2 or uppt-TH2-3) show further improvements. While the upn-EO2 is not Stokes stable, which can cause fundamental problems near steep pressure gradients (see e.g. [42]), the uppn formulations remains sta- ble. Unfortunately a direct transfer to the lTPM resp. the TPM is not possible, due to the equation structure. Replacing the total pressure by the solid pressure, enables a trans- fer of the formulation to the lTPM, but prevents the usage of equal order approximations for displacement and fluid- pressure for some parameter combinations. Further research on stable approximations using equal approximation orders Table 8 Mesh series and the resulting number of degrees of freedom (in brackets: nS ∈ DP). The meshes are characterized by the number of elements in y-direction (nelmt_y) nelmt_y 10 20 40 80 160 320 640 upn-EO2 211 411 1338 3705 13646 48081 187678 upn-mini-0 126 246 849 2412 9141 32676 128709 upn-mini-1 128 (166) 248 (326) 812 (1169) 2256 (3372) 8348 (12981) 29488 (46756) 115292 (185029) upn-TH1 170 (208) 330 (408) 1056 (1413) 2902 (4018) 10600 (15233) 37190 (54458) 144776 (214513) upn-TH2 374 (431) 734 (851) 2504 (3059) 7074 (8827) 26624 (33971) 94834 (122331) 372704 (484019) uppn-TH1-2 233 (271) 453 (531) 1461 (1818) 4029 (5145) 14773 (19406) 51933 (69201) 202421 (272158) uppn-TH2-3 498 (555) 978 (1095) 9484 (3096) 9484 (11237) 35763 (43110) 127508 (155005) 501411 (612726) uvpn-1 126 246 849 2412 9141 32676 128709 uvpn-2 209 409 1375 3861 14439 51269 201095 123 Computational Mechanics (2025) 75:1191–1212 1207 Fig. 16 Convergence over solution time per Newton iteration of displacement ũS, pressure p̃, volume-fraction nS and outflow (numerrMO) for different discretisation schemes for displacement and pressure is therefore required. A pos- sible outway might be offered by the usage of Fortin-Soulie elements [68] for the displacement. Beside the field quantities, an accurate outflow is required in the determination of the mass error. The convergence of the error of the outflow (72) numerrMO = |M̃out − M̃out,ext| M̃out,ext (79) is depicted within Fig. 16. The poor performance of approxi- mations using first order pressure approximations is obvious. Improvements are possible, either by higher order pressure approximations or by using formulations resolving the seep- age velocity. The second approach increases the accuracy of the outflow, while leading to an insufficient approximation of the pressure. Further, the second order elements for the displacement do not show the expected convergence ratio of two, as it was reported for manufactured solutions in literature [62]. However, flux based formulations are attractive, especially for the modelling of convective transport processes of dis- solved substrates. The numerical solution of the underlying advection–diffusion-reaction equation, typically requires an H(div) conforming velocity field. Higher order schemes (e.g. 123 1208 Computational Mechanics (2025) 75:1191–1212 [63], [64]) might resolve this issue, but are, due to their non- conforming nature, not always implementable. Up to this point only the spacial aspects of the discretisa- tion were analysed.While the accuracy of low order schemes (upper part of Table 7) is mainly restricted by the spa- cial resolution, higher order approximations are limited by the time discretisation, in the considered case the common implicit Euler method is used. Higher order methods would be desirable and seem to offer potential, in the construction of accurate and efficient solution procedures but are left as future work. 5 Conclusion Phase transition processes in porous materials undergoing consolidation are analysed.Homgenisationbased approaches, either following the TPM in a fully-nonlinear respectively linearised version or a phenomenological model based on Biot are utilised. Particular emphasis of this contribution is the influences of different kinematics on the time-dependent behaviour of the solution as well as its mass conserva- tivity. Therefore, a classical Terzaghi problem is modified by including phase-exchange between solid and fluid. The phases are assumed to have different densities. In order to determine the interaction of the different model parameters, non-dimensional forms of the underlying equation systems are developed. Based on an analytic solution of the problem, using the Biot equations, an inherent mass error is shown. The follow- ing findings are relevant: • Mass error and maximal displacement coincide • Themass error increases with increasing phase transition and decreasing compressibility • The mass error for cases with outflow is dominated by errors resulting from the wrong flux approximation These findings also hold true for the lTPM. Using the fully non-linear TPM, the results are mass conservative, if an appropriate discretisation scheme is used. Beside mass con- servation, models with simplified kinematics further affect the solution: • The simplified models underestimates the maximum dis- placement while overestimating the fluid pressure within the system. • Independent of the actualmass loss, simplifiedmodels for poro-elasticity with phase transition suffer in predicting the consolidation process during and followed by phase transition. This becomes especially relevant when con- solidation is neither very fast nor very slow. • Different than e.g. in classic solid mechanics there seem to be no clear limit with respect to maximal strains up to which small deformation models can be applied. It depends, beside the displacement, on the consolidation process. • The lTPM(R) model shows good approximations of the maximum displacements (compared to the TPM) espe- cially when larger displacements are present. At the same time the fluid pressure is highly overestimated and the temporal evolution of the consolidation is to fast. • The lTPM(B) model tends to underestimate the max- imum displacements compared to the TPM. Just like the other simplified model the fluid pressure is over- estimated, but not as considerably as for the lTPM(R). Compared to the other simplified models, the evolution in time is most similar to the results computed based on the TPM. For practictise relevant applications, appropriate numer- ical solution procedures are required. Comparing different numerical schemes, based on different primal variables and FE-spaces shows the following results: • Discontinuous approximations for the volume fraction increase the governing equation system and therefore the storage requirement but typically lead to faster solution procedures. • Formulations based on equal approximations orders of displacement and pressure outperform the classi- cal Taylor-Hood based upn-formulation. Higher order approximations based on the uppn-formulation show a superior efficiency. • If classical upn-formulations based on Taylor-Hood ele- ments have to be used, as e.g. the uppn-formulation is not transferable to the (l)TPM, the displacement order should be at least three. • Especially the accuracy of the higher order schemes is limited by the first order time disctetisation. Future research will cover the transfer of the here gained results onto practise relevant cases. This requires the incor- poration of temperature resp. the energy equation into the model. The idea of non-dimensionalisation in order to con- struct self-similar solutions will be applied. As classical applications in e.g. soil mechanics often deal with partially saturated porous media, the incorporation of an additional gas-phase seems reasonable. On the numerical side local mass conservation will become relevant. Within practical applications phase transition occurs locally and affects the material parameters, which will followingly become dis- continous. Classical discretisation schemes, which are not locally conservative, produce inaccurate solution in such situations, even for cases without phase transition. Higher 123 Computational Mechanics (2025) 75:1191–1212 1209 order schemes based on displacement, pressure, flux and vol- ume fraction or enrichedGalerkinmethods for displacement, pressure and volume-fraction seems to be promising. To fully utilise the potential of this methods, higher-order alternatives for the classical BDF1 method in time are required. Appendix A weak formulations Within Sect. 2.4 two alternative discretisations for poro- elastic problems were highlighted. As literature contains weak-forms only for the dimensional case, as well as for problems without phase transition, required modifications are presented in this appendix. As within Sect. 2.4, this is done for the fully non-linear TPM. Excluded from this is the uppn fomrulation, as is can be be only used for the Biot case. The uppn-formulation Choosing displacement, fluid pressure, total pressure and volume-faction as primal unknowns, requires (ũS, p̃, p̃t, nS) ∈ ( H1 )d × H1 × H1 × L2 such that a1([ũS, p̃t], vu) + a2([ũS, nS], vu) = l1(vu), b1(φ, vp) + b2(p̃, vp) = l̂2(vp), d([ũS, p̃, p̃t], vpt) = 0, c([φ, φn, nS], vn) = l̂3(vn), (A1) with φ = ( p̃ − p̃t ) / 1λ̃ S holds for all test functions (vu, vp, vpt, vn) ∈ ( H1 )d × H1 × H1 × L2. The modified right-hand-sides l̂ i , required due to the time discretisation, are defined by l̂2(vp) = b1(φn, vp) + l2(vp) , (A2) l̂3(vn) = ∫ �0 ( nS )n · vn dV , (A3) where φn and ( nS )n denote solutions of the last converged time-step.As alreadymentioned, the total-pressure based for- mulation is only applicable to the Biot problem (42). Based on the definition of the total pressure (52) the stress definition can be rewritten: ( T̃S + T̃F ) ∣∣∣∣ uppn = 2μ̃S ε̃S − p̃t I . (A4) Let now�0 be the reference configuration and � ⊆ ∂�0, the forms ai , bi , c, d and li are defined by a1 = ∫ �0 ( T̃S + T̃F ) ∣∣∣∣ uppn · ∇X̃ vu dV , a2 = ∫ �0 − [ 5 w̃FS ∣∣ B − s1 ] · vu dV , l1 = ∫ � ( T̃S + T̃F ) ∣∣∣∣ uppn · N · vu dA , (A5) b1 = ∫ �0 φ · vp dV , b2 = ∫ �0 − t̃ nFw̃FS ∣∣ B · ∇X̃ vp dV , l2 = ∫ �0 s2 · vp dV − ∫ � t̃ nFw̃FS ∣∣ B · N · vp dA , (A6) c = ∫ �0 [ nS [ φ + 1] − t̃ s3 ] · vn dV , (A7) d = ∫ �0 [ p̃ − p̃t ˜ 1λ − ∇X̃ · ũS ] · vpt dV , (A8) using φ = φ − φn . The uvpn-formulation Choosing displacement, seepage velocity, pressure and volume-faction as primal unknowns, requires (ũS, p̃, w̃, nS) ∈( H1 )d × H(div) × L2 × L2 such that a1([ũS, p̃], vu) + a2([ũS, nS], vu) = l1(vu), b1(ũS, vp) + b2(w̃, vp) = l̂2(vp), d1([ũS, w], vw) + d2(p̃, vw) = l4(vw), c1([ũS, ũnS, nS], vn) + c2(nS, vn) = l̂3(vn), (A9) holds for all test functions (vu, vw, vp, vn) ∈ ( H1 )d × H(div)×L2×L2. The modified right-hand-sides l̂ i , required due to the time discretisation, are defined by l̂2(vp) = ∫ �0 ∇X̃ · ũnS · vp dV , (A10) l̂3(vn) = c2 (( nS )n , vn ) , (A11) where ũnS and ( nS )n denote solutions of the last converged time-step. As already noted, the momentum source due to phase transition, which is typically neglected in literature, is here also neglected as is can not incorporated in this dis- cretisation. Let now �0 be the reference configuration and � ⊆ ∂�0, the forms ai , bi , di and li are defined by a1 = ∫ �0 [ P̃S + P̃F ] · ∇X̃ vu dV , a2 = ∫ �0 +S1 · vu dV , l1 = ∫ � [ P̃S + P̃F ] · N · vu dA , (A12) 123 1210 Computational Mechanics (2025) 75:1191–1212 b1 = ∫ �0 [∇X̃ · ũS − S2 ] · vp dV , b2 = ∫ �0 − t̃ ∇X̃ · w̃ · vp dV , (A13) d1 = ∫ �0 1 J̃Sk̃D [ C̃S · w̃ ] · vw dV , d2 = ∫ �0 p̃ · ∇X̃ · vw dV , l4 = ∫ � p̃ N · vw dA . (A14) The forms ci follow Eq. (51). Appendix B growth of an undrained square As it was described in Sect. 3, the current configuration of the problem stated in Fig. 3 is given by Eq. (53). Within this section solutions for displacement and volume-fraction, describing phase-transition by the fully-nonlinear TPM, are derived. The same methodology can be applied to the sim- plified models after Biot or linearisations of the TPM. For the problem at hand the non-dimensional displace- ment gradient, respectively its material time derivative are given by F̃S = [ 1 + ã 0 0 1 + ã ] , (B15) and 〈 F̃S 〉′ S = ∇X̃ 〈ũS〉′S = [ ã′ 0 0 ã′ ] . (B16) As the pressure is constant within the domain the the seepage velocity vanishes. Assuming constant material parameters, the overall mass balance (37)2 reduces to ∇x̃ · 〈ũ〉′S = J̃S c . (B17) Using Eqs. (B15) and (B16) the determinat of the deforma- tion gradient yields J̃S = [ 1 + ã ]2 , (B18) while the current gradient of the solid velocity can be expressed by ∇x̃ · 〈ũ〉′S = 2ã′ [ 1 + ã ] . (B19) Insertion into the overall mass balance results in an ODE for the non-dimensional elongation 2ã′ (1 + ã) = (1 + ã)2 c . (B20) The solution can be gained by spearing ã from its time deriva- tive ã′. The results using a(t̃) = 0 is given in Table 3. By insertion of the reduced, overallmass balance, the solid mass balance, which determines the volume fraction, simpli- fies to 〈 nS 〉′ S + nS c − 4 5 2 = 0 . (B21) The solution of this ODE, the volume fraction is spatially constant, is given by Eq. (56). Author Contributions Conceptualization, M.B., M.S., S.M.S, and TR; Methodology: M.B., M.S. and S.M.S; Implementation: M.B.; Analysis of the results: M.B., M.S. and S.M.S; Supervision: TR. and S.M.S.; Writingoriginal draft: M.B., M.S. and S.M.S; Writingreview and edit- ing, All authors. All authors have read and agreed to the published version of the manuscript. Funding Open Access funding enabled and organized by Projekt DEAL. TR thanks the Deutsche Forschungsgemeinschaft (DFG, Ger- man Research Foundation) for support via the following projects: Project-ID: 312860381 (Priority Program SPP 1886, Subproject 12); Project-ID: 390740016 (Germany’s Excellence Strategy EXC 2075/1, Project PN2-2(II)); Project-ID: 436883643 (Research Unit Programme FOR 5151 (QuaLiPerF, Project P7 Liver Lobule); Project-ID: 465194077 (Priority Programme SPP 2311, Project SimLivA); Project- ID: 327154368 (SFB 1313 Project C03 Vertebroplasty); Project-ID: 463296570 (Priority Programme SPP 1158, Antarctica), 504766766 (Project Hybrid MOR), Project-ID: 498601949 (TRR 364, Syntrac, Project B06). TR and MS are further supported by the Federal Min- istry of Education and Research (BMBF, Germany) within ATLAS by grant number 031L0304A. Data availibility statement The data that support the findings of this work are archived on the data repository od the University of Stuttgart (DaRUS) (https://doi.org/10.18419/darus-4460). All data are published under GPL 3.0 licence. For any further inquiries regarding the data or usage, please contact the corresponding author. Declarations Conflict of interest The authors declare that the researchwas conducted in the absence of any commercial or financial relationships that could be construed as a potential Conflict of interest. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adap- tation, 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 indi- cate 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, youwill need to obtain permission directly from the copy- right holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. 123 https://doi.org/10.18419/darus-4460 http://creativecommons.org/licenses/by/4.0/ http://creativecommons.org/licenses/by/4.0/ Computational Mechanics (2025) 75:1191–1212 1211 References 1. Sun W, Chen Q, Ostien JT (2014) Modeling the hydro-mechanical responses of strip and circular punch loadings on water-saturated collapsible geomaterials. Acta Geotech 9:903–934. https://doi.org/ 10.1007/s11440-013-0276-x 2. Ricken T, de Boer R (2003) Multiphase flow in a capillary porous medium. Comput Mater Sci 28(3):704–713. https://doi.org/10. 1016/j.commatsci.2003.08.032 3. Soltani K, Seyedpour SM, Ricken T, Rezazadeh G (2024) Tran- sient high-frequency sphericalwave propagation in porousmedium using fractional calculus technique. ActaMech 235(4):1845–1863. https://doi.org/10.1007/s00707-023-03780-3 4. Ricken T, Thom A, Pathak R, Thoms S, Kutschan B (2023) A continuum multi-phase and multi-component description of the freezing process of seawater to porous sea ice. https://doi.org/10. 21203/rs.3.rs-3556971/v1 5. Seyedpour S, Thom A, Ricken T (2023) Simulation of contam- inant transport through the vadose zone: A continuum mechan- ical approach within the framework of the extended theory of porous media (etpm). Water 15(2):343. https://doi.org/10.3390/ w15020343 6. Seyedpour S, Henning C, Kirmizakis P, Herbrandt S, Ickstadt K, Doherty R, Ricken T (2022) Uncertainty with varying sub- surface permeabilities reduced using coupled random field and extended theory of porous media contaminant transport models. Water 15(1):159. https://doi.org/10.3390/w15010159 7. Ateshian GA, Ricken T (2010) Multigenerational interstitial growth of biological tissues. Biomech Model Mechanobiol 9:689– 702. https://doi.org/10.1007/s10237-010-0205-y 8. Seyedpour SM, Nabati M, Lambers L, Nafisi S, Tautenhahn H- M, Sack I, Reichenbach JR, Ricken T (2021) Application of magnetic resonance imaging in liver biomechanics: a systematic review. Front Physiol 12:733393. https://doi.org/10.3389/fphys. 2021.733393 9. Seyedpour SM, Nafisi S, Nabati M, Pierce DM, Reichen- bach JR, Ricken T (2022) Magnetic resonance imaging-based biomechanical simulation of cartilage: A systematic review. J Mech Behav Biomed Mater 126:104963. https://doi.org/10.1016/ j.jmbbm.2021.104963 10. SuditschM,RickenT,WagnerA (2023) Patient-specific simulation of brain tumour growth and regression. PAMM 23(1):202200213. https://doi.org/10.1002/pamm.202200213 11. Bowen RM (1969) The thermochemistry of a reacting mixture of elastic materials with diffusion. Arch Ration Mech Anal 34:97– 127. https://doi.org/10.1007/BF00247461 12. Passman S (1977) Mixtures of granular materials. Int J Eng Sci 15(2):117–129. https://doi.org/10.1016/0020-7225(77)90027-1 13. Drumheller DS (1978) The theoretical treatment of a porous solid using a mixture theory. Int J Solid Struct 14(6):441–456. https:// doi.org/10.1016/0020-7683(78)90009-4 14. de Boer R (2012) Theory of porous media: highlights in historical development and current state, Spronger, Berlin. https://doi.org/ 10.1007/978-3-642-59637-7 15. Ehlers W, Bluhm J (2002) Porous media: theory experiments and numerical applications. Springer, Berlin. https://doi.org/10.1007/ 978-3-662-04999-0 16. HassanizadehM, GrayWG (1979) General conservation equations formulti-phase systems: 1. averagingprocedure.AdvWaterResour 2:131–144. https://doi.org/10.1016/0309-1708(79)90025-3 17. HassanizadehM, GrayWG (1979) General conservation equations for multi-phase systems: 2. mass, momenta, energy, and entropy equations. Adv Water Resour 2:191–203. https://doi.org/10.1016/ 0309-1708(79)90035-6 18. Lewis RW, Schrefler BA (1999) The finite element method in the static and dynamic deformation and consolidation of porousmedia. John Wiley & Sons, Chinchester 19. Biot MA (1941) General theory of three-dimensional consol- idation. J Appl Phys 12(2):155–164. https://doi.org/10.1063/1. 1712886 20. Schiffmann R (1971) A thermoelastic theory of consolidation. Environmental and geophysical heat transfer 21. McTigue DF (1986) Thermoelastic response of fluid-saturated porous rock. J Geophys Res: Solid Earth 91(B9):9533–9542. https://doi.org/10.1029/JB091iB09p09533 22. Pao WK, Lewis RW, Masters I (2001) A fully coupled hydro- thermo-poro-mechanical model for black oil reservoir simulation. Int J Numeric Anal Method Geomech 25(12):1229–1256. https:// doi.org/10.1002/nag.174 23. Vilchevskaya EN, Ivanova EA, Altenbach H (2014) Description of liquid–gas phase transition in the frame of continuum mechanics. Continuum Mech Thermodyn 26:221–245 24. Ehlers W, Häberle K (2016) Interfacial mass transfer during gas- liquid phase change in deformable porous media with heat transfer. Transp Porous Media 114(2):525–556. https://doi.org/10.1007/ s11242-016-0674-2 25. Ricken T, Bluhm J (2010) Modeling fluid saturated porous media under frost attack. GAMM-Mitteilungen 33(1):40–56. https://doi. org/10.1002/gamm.201010004 26. Bluhm J, Ricken T, Bloßfeld M (2011) Ice formation in porous media. Adv Ext Multifield Theor Contin 59:153–174. https://doi. org/10.1007/978-3-642-22738-7_8 27. Arasteh-Khoshbin O, Seyedpour S, Brodbeck M, Lambers L, Ricken T (2023) On effects of freezing and thawing cycles of concrete containing nano-sio 2: experimental study of material properties and crack simulation. Sci Rep 13(1):22278. https://doi. org/10.1038/s41598-023-48211-4 28. Ricken T, Schwarz A, Bluhm J (2006) A triphasic theory for growth in biological tissue-basics and applications. Materialwis- senschaft undWerkstofftechnik 37(6):446–456. https://doi.org/10. 1002/mawe.200600018 29. Ricken T, Bluhm J (2009) Evolutional growth and remodeling in multiphase living tissue. ComputMater Sci 45(3):806–811. https:// doi.org/10.1007/s00419-009-0383-1 30. Ricken T, Bluhm J (2010) Remodeling and growth of living tissue: a multiphase theory. Arch Appl Mech 80(5):453–465. https://doi. org/10.1007/s00419-009-0383-1 31. Wang H (2000) Theory of linear poroelasticity with applications to geomechanics and hydrogeology. In: Princeton University Press, Princeton. https://doi.org/10.1515/9781400885688 32. Geng Y, Yu H, McDowell G (2013) Discrete element modelling of cavity expansion and pressuremeter test. Geomech Geoeng 8(3):179–190. https://doi.org/10.1080/17486025.2012.735375 33. Griffiths D (1989) Computation of collapse loads in geomechanics by finite elements. Ingenieur-Archiv 59(3):237–244. https://doi. org/10.1007/BF00532253 34. EhlersW,Graf T, AmmannM (2004)Deformation and localization analysis of partially saturated soil. Comput Method Appl Mech Eng 193(27–29):2885–2910. https://doi.org/10.1016/j.cma.2003. 09.026 35. Lambers L, Waschinsky N, Schleicher J, König M, Taut- enhahn HM, Albadry M, Ricken T (2024) Quantifying fat zonation in liver lobules: an integrated multiscale in silico model combining disturbed microperfusion and fat metabolism via a continuum biomechanical bi-scale, tri-phasic approach. Biomech Model Mechanobiol 23(2):631–653. https://doi.org/10. 1007/s10237-023-01797-0 36. Sanei M, Duran O, Devloo PR, Santos ES (2022) Evaluation of the impact of strain-dependent permeability on reservoir produc- tivity using iterative coupled reservoir geomechanical modeling. 123 https://doi.org/10.1007/s11440-013-0276-x https://doi.org/10.1007/s11440-013-0276-x https://doi.org/10.1016/j.commatsci.2003.08.032 https://doi.org/10.1016/j.commatsci.2003.08.032 https://doi.org/10.1007/s00707-023-03780-3 https://doi.org/10.21203/rs.3.rs-3556971/v1 https://doi.org/10.21203/rs.3.rs-3556971/v1 https://doi.org/10.3390/w15020343 https://doi.org/10.3390/w15020343 https://doi.org/10.3390/w15010159 https://doi.org/10.1007/s10237-010-0205-y https://doi.org/10.3389/fphys.2021.733393 https://doi.org/10.3389/fphys.2021.733393 https://doi.org/10.1016/j.jmbbm.2021.104963 https://doi.org/10.1016/j.jmbbm.2021.104963 https://doi.org/10.1002/pamm.202200213 https://doi.org/10.1007/BF00247461 https://doi.org/10.1016/0020-7225(77)90027-1 https://doi.org/10.1016/0020-7683(78)90009-4 https://doi.org/10.1016/0020-7683(78)90009-4 https://doi.org/10.1007/978-3-642-59637-7 https://doi.org/10.1007/978-3-642-59637-7 https://doi.org/10.1007/978-3-662-04999-0 https://doi.org/10.1007/978-3-662-04999-0 https://doi.org/10.1016/0309-1708(79)90025-3 https://doi.org/10.1016/0309-1708(79)90035-6 https://doi.org/10.1016/0309-1708(79)90035-6 https://doi.org/10.1063/1.1712886 https://doi.org/10.1063/1.1712886 https://doi.org/10.1029/JB091iB09p09533 https://doi.org/10.1002/nag.174 https://doi.org/10.1002/nag.174 https://doi.org/10.1007/s11242-016-0674-2 https://doi.org/10.1007/s11242-016-0674-2 https://doi.org/10.1002/gamm.201010004 https://doi.org/10.1002/gamm.201010004 https://doi.org/10.1007/978-3-642-22738-7_8 https://doi.org/10.1007/978-3-642-22738-7_8 https://doi.org/10.1038/s41598-023-48211-4 https://doi.org/10.1038/s41598-023-48211-4 https://doi.org/10.1002/mawe.200600018 https://doi.org/10.1002/mawe.200600018 https://doi.org/10.1007/s00419-009-0383-1 https://doi.org/10.1007/s00419-009-0383-1 https://doi.org/10.1007/s00419-009-0383-1 https://doi.org/10.1007/s00419-009-0383-1 https://doi.org/10.1515/9781400885688 https://doi.org/10.1080/17486025.2012.735375 https://doi.org/10.1007/BF00532253 https://doi.org/10.1007/BF00532253 https://doi.org/10.1016/j.cma.2003.09.026 https://doi.org/10.1016/j.cma.2003.09.026 https://doi.org/10.1007/s10237-023-01797-0 https://doi.org/10.1007/s10237-023-01797-0 1212 Computational Mechanics (2025) 75:1191–1212 Geomech Geophys Geo-Energy Geo-Resour 8(2):54. https://doi. org/10.1007/s40948-022-00344-y 37. Barbeiro S (2019) Poroelasticity with deformation dependent per- meability. Progress in industrial mathematics at ECMI 2018. Springer, Berlin, pp 389–395 38. Pierce DM, Ricken T, Holzapfel GA (2013) A hyperelastic biphasic fibre-reinforced model of articular cartilage considering distributed collagen fibre orientations: continuum basis, computa- tional aspects and applications. Comput Method Biomech Biomed Eng 16(12):1344–1361. https://doi.org/10.1080/10255842.2012. 670854 39. Murad MA, Loula AFD (1992) Improved accuracy in finite ele- ment analysis of Biot’s consolidation problem. Comput Method Appl Mech Eng 95(3):359–382. https://doi.org/10.1016/0045- 7825(92)90193-N 40. Murad MA, Loula AFD (1994) On stability and convergence of finite element approximations of Biot’s consolidation problem. Int J Numeric Method Eng 37(4):645–667. https://doi.org/10.1002/ nme.1620370407 41. Stokes IA, Chegini S, Ferguson SJ, Gardner-Morse MG, Iatridis JC, Laible JP (2010) Limitation of finite element analysis of poroe- lastic behavior of biological tissues undergoing rapid loading. Ann Biomed Eng 38:1780–1788. https://doi.org/10.1007/s10439-010- 9938-0 42. Bertrand F, Brodbeck M, Ricken T (2022) On robust discretiza- tion methods for poroelastic problems: Numerical examples and counter-examples. Ex Counterex 2:100087. https://doi.org/10. 1016/j.exco.2022.100087 43. Haga JB,OsnesH,LangtangenHP (2012)On the causes of pressure oscillations in low-permeable and low-compressible porousmedia. Int J Numeric Anal Method Geomech 36(12):1507–1522. https:// doi.org/10.1002/nag.1062 44. Berger L, Bordas R, Kay D, Tavener S (2017) A stabilized finite element method for finite-strain three-field poroelasticity. Comput Mech 60:51–68. https://doi.org/10.1007/s00466-017-1381-8 45. Liu J, Mu L, Ye X (2011) A comparative study of locally conser- vative numerical methods for darcy’s flows. Procedia Comput Sci 4:974–983. https://doi.org/10.1016/j.procs.2011.04.103 46. Scovazzi G, Wheeler MF, Mikelić A, Lee S (2017) Analytical and variational numerical methods for unstable miscible displacement flows in porous media. J Comput Phys 335:444–496. https://doi. org/10.1016/j.jcp.2017.01.021 47. Choo J, Lee S (2018) Enriched Galerkin finite elements for cou- pled poromechanicswith localmass conservation. ComputMethod ApplMechEng 341:311–332. https://doi.org/10.1016/j.cma.2018. 06.022 48. Oyarzúa R, Ruiz-Baier R (2016) Locking-free finite element meth- ods for poroelasticity. SIAM J Numeric Anal 54(5):2951–2973. https://doi.org/10.1137/15M1050082 49. Lee JJ,MardalK-A,Winther R (2017) Parameter-robust discretiza- tion and preconditioning ofBiot’s consolidationmodel. SIAMJSci Comput 39(1):1–24. https://doi.org/10.1137/15M1029473 50. Drieschner M, Matthies HG, Hoang T-V, Rosić B, Ricken T, Henning C, Ostermeyer G-P, Müller M, Brumme S, Srisupat- tarawanit T, Weinberg K, Korzeniowski TF (2019) Analysis of polymorphic data uncertainties in engineering applications. GAMM-Mitteilungen 42(2):201900010. https://doi.org/10.1002/ gamm.201900010 51. Ricken T, Schröder J, Bluhm J, Maike S, Bartel F (2022) The- oretical formulation and computational aspects of a two-scale homogenization scheme combining the tpm and fe 2 method for poro-elastic fluid-saturated porous media. Int J Solid Struct 241:111412. https://doi.org/10.1016/j.ijsolstr.2021.111412 52. Ehlers W (1996) Grundlegende Konzepte in der Theorie Poröser Medien. Technische Mechanik 16:63–76 53. Truesdell C (1984) Rational thermodynamics. Springer, NewYork. https://doi.org/10.1007/978-1-4612-5206-1 54. Simo JC, Pister KS (1984) Remarks on rate constitutive equa- tions for finite deformation problems: computational implications. Comput Method Appl Mech Eng 46(2):201–215. https://doi.org/ 10.1016/0045-7825(84)90062-8 55. Diebels S, Ehlers W, Markert B (2001) Neglect of the fluid extra stresses in volumetrically coupled solid-fluid problems. ZAMM-J Appl Math Mech/Zeitschrift für Angewandte Math- ematik und Mechanik 81(S3):521–522. https://doi.org/10.1002/ zamm.20010811540 56. Ehlers W, Morrison M, Schröder P, Stöhr D, Wagner A (2021) Multiphasic modelling and computation of metastatic lung-cancer cell proliferation and atrophy in brain tissue based on experimental data.BiomechModelMechanobiol 21(1):277–315. https://doi.org/ 10.1007/s10237-021-01535-4 57. Terzaghi K (1943) Theoretical soil mechanics. JohnWiley & Sons, Ltd, New York. https://doi.org/10.1002/9780470172766 58. Soderberg LO (1962) Consolidation theory applied to foundation pile time effects. Géotechnique 12(3):217–225. https://doi.org/10. 1680/geot.1962.12.3.217 59. OsmanAS, RandolphMF (2010) Response of a solid infinite cylin- der embedded in a poroelastic medium and subjected to a lateral load. Int J Solid Struct 47(18):2414–2424. https://doi.org/10.1016/ j.ijsolstr.2010.05.001 60. Ricken T, Schwarz A, Bluhm J (2007) A triphasic model of trans- versely isotropic biological tissue with applications to stress and biologically induced growth. Comput Mater Sci 39(1):124–136. https://doi.org/10.1016/j.commatsci.2006.03.025 61. Ehlers W, Markert B, Röhrle O (2009) Computational continuum biomechanics with application to swelling media and growth phe- nomena. GAMM-Mitteilungen 32(2):135–156. https://doi.org/10. 1002/gamm.200910013 62. Mardal K-A, Rognes ME, Thompson TB (2021) Accurate dis- cretization of poroelasticity without darcy stability. BIT Numeric Math. https://doi.org/10.1007/s10543-021-00849-0 63. Hong Q, Kraus J (2018) Parameter-robust stability of classi- cal three-field formulation of Biot’s consolidation model. ETNA - Electron Trans Numeric Anal 48:202–226. https://doi.org/10. 1553/etna_vol48s202 64. Kraus J, Lederer PL, Lymbery M, Schöberl J (2021) Uni- formlywell-posedhybridizeddiscontinuous galerkin/hybridmixed discretizations for biot’s consolidation model. Comput Method Appl Mech Eng 384:113991. https://doi.org/10.1016/j.cma.2021. 113991 65. Stickle MM, Pastor M (2018) A practical analytical solution for one-dimensional consolidation. Géotechnique 68(9):786–793. https://doi.org/10.1680/jgeot.16.P.268 66. Baratta IA, Dean JP, Dokken JS, Habera M, Hale JS, Richard- son CN, Rognes ME, Scroggs MW, Sime N, Wells GN (2023) DOLFINx: The next generation FEniCS problem solving environ- ment. Zenodo. https://doi.org/10.5281/zenodo.10447666 67. Amestoy PR, Duff IS, L’Excellent JY, Koster J (2001) A fully asyn- chronousmultifrontal solver using distributed dynamic scheduling. SIAM J Matrix Anal Appl 23(1):15–41 68. Fortin M, Soulie M (1983) A non-conforming piecewise quadratic finite element on triangles. Int J Numeric Method Eng 19(4):505– 520. https://doi.org/10.1002/nme.1620190405 Publisher’s Note Springer Nature remains neutral with regard to juris- dictional claims in published maps and institutional affiliations. 123 https://doi.org/10.1007/s40948-022-00344-y https://doi.org/10.1007/s40948-022-00344-y https://doi.org/10.1080/10255842.2012.670854 https://doi.org/10.1080/10255842.2012.670854 https://doi.org/10.1016/0045-7825(92)90193-N https://doi.org/10.1016/0045-7825(92)90193-N https://doi.org/10.1002/nme.1620370407 https://doi.org/10.1002/nme.1620370407 https://doi.org/10.1007/s10439-010-9938-0 https://doi.org/10.1007/s10439-010-9938-0 https://doi.org/10.1016/j.exco.2022.100087 https://doi.org/10.1016/j.exco.2022.100087 https://doi.org/10.1002/nag.1062 https://doi.org/10.1002/nag.1062 https://doi.org/10.1007/s00466-017-1381-8 https://doi.org/10.1016/j.procs.2011.04.103 https://doi.org/10.1016/j.jcp.2017.01.021 https://doi.org/10.1016/j.jcp.2017.01.021 https://doi.org/10.1016/j.cma.2018.06.022 https://doi.org/10.1016/j.cma.2018.06.022 https://doi.org/10.1137/15M1050082 https://doi.org/10.1137/15M1029473 https://doi.org/10.1002/gamm.201900010 https://doi.org/10.1002/gamm.201900010 https://doi.org/10.1016/j.ijsolstr.2021.111412 https://doi.org/10.1007/978-1-4612-5206-1 https://doi.org/10.1016/0045-7825(84)90062-8 https://doi.org/10.1016/0045-7825(84)90062-8 https://doi.org/10.1002/zamm.20010811540 https://doi.org/10.1002/zamm.20010811540 https://doi.org/10.1007/s10237-021-01535-4 https://doi.org/10.1007/s10237-021-01535-4 https://doi.org/10.1002/9780470172766 https://doi.org/10.1680/geot.1962.12.3.217 https://doi.org/10.1680/geot.1962.12.3.217 https://doi.org/10.1016/j.ijsolstr.2010.05.001 https://doi.org/10.1016/j.ijsolstr.2010.05.001 https://doi.org/10.1016/j.commatsci.2006.03.025 https://doi.org/10.1002/gamm.200910013 https://doi.org/10.1002/gamm.200910013 https://doi.org/10.1007/s10543-021-00849-0 https://doi.org/10.1553/etna_vol48s202 https://doi.org/10.1553/etna_vol48s202 https://doi.org/10.1016/j.cma.2021.113991 https://doi.org/10.1016/j.cma.2021.113991 https://doi.org/10.1680/jgeot.16.P.268 https://doi.org/10.5281/zenodo.10447666 https://doi.org/10.1002/nme.1620190405 Phase transition in porous materials: effects of material parameters and deformation regime on mass conservativity Abstract 1 Introduction 2 Material methods 2.1 Governing equations 2.2 Nondimensionalisation 2.3 Pull-back and linearisation 2.4 Weak forms and numerical treatment 3 A first glance on mass losses in linearised poro-elasticity 4 Results and discussion 4.1 Problem and analytic solution 4.2 Influence of model parameters on the mass loss 4.3 Influence of the model on the solution 4.4 Efficiency of numerical solution procedures 5 Conclusion Appendix A weak formulations Appendix B growth of an undrained square References