Pore Pressure Evolution and F-T Fatigue of Concrete: A Coupled THM-F Phase-Field Modeling Approach

1  Introduction

As one of the most widely used construction materials globally, concrete plays a crucial role in infrastructure development. Its durability directly affects structural safety and long-term maintenance costs. In cold regions, such as high-altitude or high-latitude areas, freeze-thaw (F-T) damage represents a major threat to the service life of concrete structures [1]. Extensive research has confirmed significant degradation in both strength and dynamic elastic modulus of concrete under F-T conditions. With increasing F-T cycles, progressive deterioration of the pore structure occurs, leading to a higher concentration of macro-pores [2]. This process triggers the development of dense microcrack networks that propagate and coalesce, eventually resulting in structural fracture [2,3]. Given the severe impact of F-T damage, understanding its underlying mechanisms is essential for extending the service life of concrete structures.

Recent experimental research in F-T damage has shifted from macroscopic performance evaluations toward microstructural characterization, with a particular emphasis on the evolution of pore structure. Early studies predominantly relied on conventional metrics such as mass loss and relative dynamic elastic modulus to assess F-T deterioration levels [2], yet these approaches offered limited insight into the intrinsic damage mechanisms. Advances in non-destructive testing techniques, including X-ray computed tomography (CT), low-field nuclear magnetic resonance (LF-NMR), and scanning electron microscopy, have enabled detailed visualization of pore structure evolution during F-T cycles. For instance, Zhou et al. [2] employed X-ray CT to document increased macro-pore content and the formation interconnected pore network under F-T cyclic conditions. Zhang et al. [4] developed a porosity-based damage model using LF-NMR to quantitative correlate microstructural changes with macroscopic mechanical degradation in concrete with varying water-cement ratios. Kristensen et al. [5] further utilized micro-scale tomography to track crack initiation and propagation within the interfacial transition zone (ITZ). These studies collectively demonstrate that mesoscale structural evolution, including changes in pore size distribution, connectivity, and microcracking, directly governs the macroscopic mechanical degradation induced by F-T cycles.

Despite these advances, experimental methods are inherently constrained by long durations, high costs, and limited operational conditions, which impede comprehensive capture of the entire damage evolution process. These limitations have motivated the development of theoretical and numerical models to elucidate damage mechanisms and improve engineering predictability.

F-T damage in concrete is a multifactorial coupled process. Theoretical models bridge macro-behavior with micro-mechanisms, notably through Powers’s hydrostatic pressure hypothesis [6] and Powers and Helmuth’s improved osmotic pressure theory [7]. Both attribute damage primarily to expansive stress induced by pore ice-water phase transition and osmotic pressure caused by unfrozen water migration. Subsequent theories, including the pore crystallization hypothesis [8] and bond detachment theory [9], further elucidate the mechanisms of ice formation and surface scaling. Setzer and Auberg [10] derived thermodynamic-based pressure equations for saturated porous media during freezing, while Penttala [11] established an analogous theoretical model and verified the expansion formula for triple-phase pore systems. Coussy’s mechanical model of porous media [12] provided a foundational framework linking microstructural evolution to macroscopic mechanical behavior.

Recognizing F-T as a coupled multiphysics process, recent research has focused on integrated modeling approaches. Within the continuum damage mechanics (CDM) framework, F-T damage is defined as the evolution of internal damage variables with cycles or time. By incorporating parameters such as temperature, saturation, and ice-water phase transition, several coupled models have been proposed. For instance, Liu et al. [13] established a comprehensive thermo-hydro-mechanical (THM) model for fractured porous media. This model fully couples the constitutive behavior of elastoplastic materials with a continuum damage mechanics framework. Another study by Niu et al. [14] proposed a damage constitutive model for concrete that incorporates the coupling effect of F-T cycles and external loads, utilizing homogenization theory. Yu et al. [15] presented a fully coupled THM numerical implementation for solid-water-ice systems employing the Material Point Method. Although CDM provides a fundamental theoretical framework for understanding concrete F-T damage mechanisms, it falls short in capturing microstructural evolution, strong physical couplings, and crack propagation, particularly under F-T conditions where cumulative micro-damage is critical for durability assessment.

Numerical methods have demonstrated strong capabilities in simulating F-T damage, with high-resolution meshing improving accuracy for localized cracking and nonlinear stress fields. The team of Professor Ueda and Wang [16] integrated advanced constitutive relations into the Rigid Body Spring Model, effectively capturing the F-T damage evolution and providing insight into crack propagation and multiscale coupling. Jun et al. [17] formulated a Discrete Element Method incorporating pore-water phase transition and ice-crystal expansion. Wu et al. [18] and Li et al. [19] applied Peridynamics modeling for F-T failure process. Zheng et al. [20] established an F-T damage model using the Material Point Method considering pore-water phase change.

Among numerical techniques, the Finite Element Method (FEM) remains prominent due to its stability and versatility. Olsen [21] developed a 2-D FEM model coupling heat and moisture transport to simulate temperature, ice content, and pore pressure in saturated concrete. Building upon this framework, Zuber and Marchand [22] further incorporated mechanical responses, thereby enhancing the model’s capability to capture freeze-induced expansion. Subsequently, based on Coussy’s poromechanics theory, Duan proposed a comprehensive THM coupling model [23] that accurately simulates temperature fields, pore pressure, and stress evolution.

However, concrete’s inherent heterogeneity (aggregate, mortar, ITZ) [18,24] leads to high stochasticity in crack propagation, requiring high simulation accuracy. High-resolution meshing required for crack initiation modeling escalates computational cost with increasing F-T cycles, and the mesh sensitivity can cause numerical instability. Current approaches fail to fully capture the micro-damage accumulation and rely heavily on empirical calibration, limiting their predictive capability for long-term performance and durability assessment.

The phase-field method (PFM) offers a promising alternative for simulating F-T damage. It represents fracture through a smooth phase-field variable that eliminates the need for explicit crack tracking or mesh updating. The bandwidth of the damage zone controlled by an intrinsic length-scale parameter. PFM can capture the initiation and propagation of cracks in heterogeneous concrete with reduced mesh sensitivity and improved numerical stability. Its inherent compatibility with multiphysics modeling provides a solid basis for simulating coupled F-T damage mechanisms [25].

Under F-T conditions, concrete damage driven by hydraulic or osmotic pressure can be regarded as a three-dimensional fatigue process under cyclic tensile stress [26]. Repeated freezing and thawing induce periodic stress fluctuations due to water-ice phase change, leading to damage accumulation at microstructural defects such as pores, microcracks, and the ITZ. Progressive micro-damage eventually coalesces into macrocracks and material degradation [26].

Fatigue modeling has evolved from empirical approaches (e.g., S-N curves [27]) toward mechanical theories. Palmgren and Miner’s linear cummulative damage rule [28,29] and Griffith’s energy-based fracture theory [30] laid early foundation. Irwin’s stress intensity factor [31,32] and Paris’ law [30] enabled quantitative prediction of fatigue crack growth rates. Recently, variational phase-field fracture model [33] has substantially advanced fatigue simulation. Based on Griffith’s energy conservation principle [30], Francfort and Marigo proposed a variational formulation for crack evolution [34], numerically implemented by Bourdin et al. using the PFM [35]. This framework was extension by Lo et al. [36], who incorporated viscous dissipation and a modified J-integral to reproduce the behavior of Paris law, while Alessi et al. [37] introduced history-dependent dissipation for cyclic loading. Carrara et al. [38] further enhanced the model’s capacity to predict the fatigue crack paths through load-history-dependent degradation function. These advances reduce reliance on empirical parameters and enhance physical consistency and predictive accuracy.

Building on advancements in phase-field methods and multiphysics modeling for concrete F-T damage, this study develops a strongly coupled thermo-hydro-mechanical fatigue (THM-F) model. While previous models have successfully integrated thermal, hydraulic, and mechanical fields, and some have incorporated phase-field approaches for fracture simulation, a key limitation remains the treatment of the material’s pore structure as static or its evolution as empirically prescribed. This study distinguishes itself by explicitly integrating physics-based pore structure evolution into the coupled THM framework. Specifically, the model incorporates a mechanistic description of how cumulative micro-damage progressively alters the pore size distribution and permeability over F-T cycles. By dynamically linking the evolving pore structure, characterized by the PSD curve, to the permeability coefficient, this work enables a more realistic simulation of the F-T damage interaction cycle-by-cycle. This integration of real-time microstructural evolution within a full THM-F framework represents a significant step forward in achieving high-fidelity prediction of F-T damage.

The remainder of this paper is organized as follows. Section 2 outlines the variational phase-field framework for fatigue. Section 3 presents the coupled multi-physics model for F-T processes, including pore structure deterioration, permeability increase, pore pressure computation, the coupling between phase-field and the THM system. Section 4 details the numerical implementation using finite element discretization and time-integration algorithms. Section 5 provides validation against experimental data, parameter determination, and discussion of damage mechanisms, microstructure degradation, hysteresis in ice content, effects of water-cement ratio, and air-void size distribution analyses for frost resistance. Finally, Section 6 summarizes the principal findings and conclusions.

2  Phase-Field Variational Framework for Fatigue

Building upon the variational theory established by Francfort and Marigo [34], which extends Griffith’s theory, the total energy functional consists of three primary components:

Ψ(u,ϕ)=∫Ωψd(ϕ)dΩ⏟Ψd+∫ΩψedΩ⏟Ψe−∫ΩαBpef(∇⋅u)dΩ−∫Ωb⋅udΩ−∫∂Ωtf¯⋅udΩ⏟Ψexternal(1)

where Ψd represents the dissipation energy, with ψd as the energy density related to dissipation, Ψe represents the elastic energy, with ψe as the energy density corresponding to elasticity; and Ψexternal accounts for the external work, including the effects of effective pore pressure pef, body forces b, and boundary forces f¯.

The damage evolution in concrete under F-T cycles originates from cyclic tensile stress accumulation induced by pore water phase transition to ice. Based on energy dissipation principles, a fatigue degradation function f(α¯(t)) is introduced to dynamically regulate the critical energy release rate [38], that is,

GcF=f(α¯(t))Gc(2)

where Gc is critical energy release rate, α¯ is a variable related to accumulated fatigue history, which can reflect the potential behavioral change under cyclic loading conditions [33]. This study adopts a mean-load-independent fatigue history variable following Carrara et al. [38], which assumes material fatigue life remains unaffected by mean load variations within a cycle. The variable evolves as:

α¯(x,t)=∫0tH(αα˙)|α˙|dτ(3)

where α is the current free energy density, that is, α=g(ϕ)ψe+(ε). H(⋅) is the Heaviside function, which ensures that damage accumulation occurs only when the effective tensile strain energy rate ψ˙e+≥0 [38]. Choose the following fatigue degradation function as

f(α¯(t))={1      if  α¯(t)≤αT[1−κlog⁡(α¯(t)αT)]2 if   αT≤α¯(t)≤αT101/(2κ)0.25 if  α¯(t)≥αT101/(2κ)(4)

where αT is the threshold parameter defined as αT=1/2Eε^2, which incorporates the characteristic strain ε^=Gc/(3l0E), synthetically characterizing fracture toughness and phase-field regularization effects [38]. The fatigue degradation function f(α¯) plays a crucial role in characterizing material toughness reduction during cyclic loading, satisfying a piecewise continuously differentiable condition. As illustrated in Fig. 1, the fatigue degradation function exhibits distinct curvature governed by the slope parameter κ∈(0,1).

Figure 1: Fatigue degradation function

Thus, the dissipation energy density ψd considering the fatigue effect can be written as:

ψd(ϕ,∇ϕ|α¯)=f(α¯)Gc(ϕ22l0+l02|∇ϕ|2)(5)

where ϕ is an order parameter reflecting material damage, ϕ=0 and ϕ=1 correspond to intact and fully fractured states, respectively. And l0 is a regularization length parameter, as shown in Fig. 2.

Figure 2: Sharp and diffusive fractures in poro-elastic medium: (a) Sharp fracture; (b) Diffusive fracture

Considering that the damage during the F-T cycles of concrete is mainly caused by tensile cracking due to frost heave, the elastic energy density ψe is formulated as

ψe(ε,ϕ)=g(ϕ)ψe+(ε)+ψe−(ε)(6)

where the degradation function g(ϕ), constrained by monotonicity (g′(ϕ)≤0) and convexity (g′′(ϕ)≥0). A widely adopted form g(ϕ)=(1−ϕ)2+kζ enables a continuous transition from intact material to complete fracture with kζ=10−9 preventing numerical singularities [39]. The strain tensor is decoupled via spectral decomposition ε±=∑i⟨εi⟩±ni⊗ni, where the bracket operator ⟨∗⟩±=(∗±|∗|)/2 selectively activates tensile components [40], ensuring crack propagation driven solely by tensile strains.

The constitutive model couples material properties through Lamé parameters, with strain energy components [40]:

ψe±(ε)=λ2⟨tr(ε)⟩±2+μtr(ε±2)(7)

where the Lamé modulus λ=Ev(1+v)(1−2v) and shear modulus μ=E2(1+v) are derived from Young’s modulus E and Poisson’s ratio v, and tr(⋅) is the trace operator.

By incorporating fracture energy and elastic energy into the total energy functional, the multi-physically coupled energy functional can be formulated as:

Ψ(u,ϕ)=∫Ω∫0tf(α¯(τ))Gc(ϕ22l0+l02|∇ϕ|2)dτdΩ⏟Ψd+∫Ω[g(ϕ)ψe+(εe)+ψe−(εe)]dΩ⏟Ψe−∫ΩαBpef(∇⋅u)dΩ−∫Ωb⋅udΩ−∫∂Ωtf¯⋅udS⏟Ψexternal(8)

The evolution of the phase-field variable is governed by the principle of energy minimization, representing the balance between fracture energy dissipation and elastic energy release. The governing equation is derived by seeking the stationarity of the total energy functional with respect to ϕ. Specifically, the first variation δΨ/δϕ=0 is calculated, considering the arbitrariness of the variational parameter δϕ. This variational process yields the following Euler-Lagrange equation [41]:

[2l0(1−kζ)Gcψe++f(α¯(t))]ϕ−l02[∇f(α¯)∇ϕ+f(α¯(t))∇2ϕ]=2l0(1−kζ)Gcψe+(9)

To ensure crack irreversibility in accordance with the second law of thermodynamics. The maximum historical value of the energy driving force is introduced [39]:

H=maxτ∈[0,t]ψe+(ε(x,τ))(10)

The complete governing equation for fatigue behavior in the fracture phase-field can be expressed as:

[2l0(1−kζ)GcH+f(α¯(t))]ϕ−l02[∇f(α¯)∇ϕ+f(α¯(t))∇2ϕ]=2l0(1−kζ)GcH(11)

With the corresponding boundary condition:

(l0f(α¯(t))Gc∇⋅ϕ)⋅n=0on ∂Ω(12)

Having established the variational phase-field framework for fatigue damage in Section 2, we now extend this foundation to incorporate multiphysics coupling essential for modeling F-T behavior in concrete. The phase-field variable ϕ, which characterizes damage evolution, will be coupled with thermal, hydraulic, and mechanical fields to capture the complex interactions during F-T cycles. This integration enables the simulation of how temperature changes drive pore water phase transitions, how resulting pressure developments induce mechanical stresses, and how these stresses ultimately lead to fatigue damage accumulation through the phase-field variable.

3  Multiphysics Coupling Framework

A critical aspect of F-T damage in concrete is the tight coupling between its thermal, hydraulic, and mechanical response mechanisms. The pore structure characteristics of saturated concrete play a decisive role in the THM coupling mechanisms during F-T cycles. Accordingly, this section first introduces the PSD curve, which characterizes the pore structure features influencing ice crystallization and phase transition. Subsequently, based on the theory of pore ice formation, the calculated model of ice saturation are derived for the saturated concrete under F-T conditions, along with models describing the evolution of porosity, permeability coefficient, and pore pressure throughout the process. Finally, these evolution models are incorporated into a THM coupling framework, leading to the establishment of a THM-F coupled finite element model designed to simulate the damage progression in concrete subjected to F-T cycles.

3.1 Pore Structure Characterization and the Porosity Evolution under F-T Condition

The predictive capability of the proposed model derives from its ability to extrapolate pore structure evolution based on the initial PSD characteristics and fundamental physical principles, rather than relying on empirically measured PSD data at specific cycle numbers. The PSD function φ(Di) is a critical parameter that characterizes the volumetric fraction of pores withdiameter Di within the porous medium. The morphological features of the PSD exert deterministic control over the ice-water phase transition behavior.

In this study, the PSD curve was quantitatively obtained through via Mercury Intrusion [42], a technique that infers the pore structure from the relationship between intrusion pressure and volume, based on the surface tension and incompressibility of mercury. As illustrated in Fig. 3, the PSD curve for the non-air-entrained concrete used in this work, validated against standardized test data from Reference [43], is mathematically represented by the following rational fitting function:

φ(x)=p1x5+p2x4+p3x3+p4x2+p5x+p6x4+q1x3+q2x2+q3x+q4(13)

where x=log10⁡Di. The fitting parameters p1−p6 and q1−q4 are listed in Table 1. These high-precision fitting results provide a reliable foundation for subsequent pore water phase transition analysis.

Figure 3: The pore size-distribution curve

According to the Gibbs-Thomson thermodynamic theory, the freezing point depression of pore water nonlinearly depends on the pore radius [11]. For saturated pores, this relationship is influenced by capillary effects described by the Young-Laplace equation, which states that the capillary pressure is inversely proportional to the pore diameter [43]. During the initial freezing stage, larger pores with greater curvature radii exhibit lower ice nucleation energy barriers, leading to preferential ice formation. As the thermal gradient intensifies, the phase transition front progressively propagates toward smaller pores, as shown in Fig. 4. The critical pore diameter Dcr serves as the key threshold parameter determining cryo-expansion susceptibility, physically demarcating the pore network into cryo-active zones (Di>Dcr), where freezing occurs and generates expansive strain, and cryo-inert zones (Di<Dcr), where water remains unfrozen due to supercooling effects [11,44,45]. Its mathematical formulation is expressed as:

Dcr=−2ΥcLT0ρwL(Tf−T0)+1.97|Tf−T0|−13(14)

where, ΥcL=39×10−3 N/m denotes the ice-water interface tension. L=333.5 kJ/kg represents the latent heat of fusion [23]. The first term characterizes the thermodynamic inhibition of ice nucleation by interface curvature effects, while the second term elucidates the exponential correlation between unfrozen water film thickness (δ) and supercooling degree (ΔT) [44]. The cryo-active zones, serving as transport pathways for ice crystal growth, synergistically regulate the microscale mechanical response mechanisms of F-T damage with the ice-suppression effects in cryo-inert zones.

Figure 4: Schematic illustration of ice crystal growth morphology within pore structures (the ice crystal tip radius Dtip, dictated by the Gibbs-Thomson relationship)

Based on the physical definition of the critical pore diameter Dcr, the cumulative volumetric fraction of pores noncryo-active zones governs the ice content at freezing temperature T. By integrating the PSD function φ(Di), the ice saturation can be mathematically expressed as:

scF=∫DcrDmaxφ(Di)dDi∫DminDmaxφ(Di)dDi=∫DcrDmaxφ(Di)dDin(15)

where, n denotes the total porosity.

Experimental studies demonstrate that geometric disparities between pore throats and bodies induce a ‘F-T hysteresis’ effect, leading to asymmetric ice content evolution between freezing and thawing, where the freezing point is controlled by throat dimensions and the melting point is governed by body dimensions. To quantify this phenomenon, the ice content during thawing is reformulated as a temperature-dependent piecewise function [46]:

scM={∫D(Tmin)Dcr(Tmin)φ(Di)dDin,Tmin<T<ϖTmin∫D(Tmin)Dcr(ϖT)φ(Di)dDin,ϖTmin<T<0(16)

where the shape factor ϖ=Tmelt/Tfreeze represents the ratio of melting point to freezing point. According to experimental studies by Sun and Scherer [47], the ϖ value for cementitious materials typically lie in the range 0.1~0.5, indicating that the melting point of pores at identical scales consistently exceeds the freezing point, as shown in Fig. 5.

Figure 5: The evolution of ice saturation during F-T cycles

The degradation of concrete under F-T cycling originates from progressive damage to its microporous structure. Experimental studies reveal a two-scale deterioration pattern: at the macroscopic level, the total porosity increases with cycle number, while microscopically, the PSD shifts toward larger diameters. To capture this dual-scale evolution, the present study develops a microstructure-informed model that links phase-change-induced deformation to pore growth.

In the model, individual pores are idealized as cylindrical inclusions—a simplification supported by microstructural studies showing that pores in cementitious materials are predominantly elongated and interconnected [48]. This approach allows applying poromechanics theory to quantify eigenstrains induced by water-ice phase change [49]. During freezing, pore water is assumed to expand uniformly in all directions. The resulting deformation is partly reversible upon thawing, but a fraction accumulates as irreversible damage, leading to pore enlargement. The change in pore size after one F-T cycle is expressed through the recursive relationship:

DiN=DiN−1{1+εr}=DiN−1{1+[1−RE(DiN−1)]αEεv∗}(17)

here, DiN is the diameter of the ith pore after N cycles; αE is an expansion amplification factor derived from Eshelby’s solution for a cylindrical inclusion, given by αE=(1+v)/[3(1−v)], where v is Poisson’s ratio of concrete; and εv∗ is the volumetric eigenstrain due to water-ice transition. The term RE(D) represents the elastic recovery coefficient, describing the fraction of expansive strain that is recovered during thawing.

Not all pores recover equally. Larger pores experience more permanent expansion due to weaker constraint from the surrounding matrix and greater interaction with neighboring pores. Experimental observations suggest that deformation recovery decreases as pore size increases [2]. To reflect this size-dependent behavior, a segmented recovery function is introduced:

Re(Di)={ log(Di)/log(Dmax)Di>100 μm0.9960 μm < Di < 100 μm​log(Di)/log(60)10 μm < Di < 60 μm0.984 μm < Di < 10 μm0.97 Dcr μm < Di < 4 μm (18)

This function assigns higher recovery to smaller pores (0.99 for pores below 4 µm) and lower recovery to larger ones (0.97 for pores 60–100 µm), reflecting the increasing irreversibility of deformation with pore size.

Based on the pore size update rule in Eq. (19), the evolution of the PSD curve after N cycles is derived as:

φN(DiN)=φ0(Di0)(DiNDi0)2∑kφ0(Dk0)(DkNDk0)2(19)

here, φ0(Di0) is the initial volume fraction of pores with diameter Di0. The square ratio accounts for the cross-sectional area change assuming cylindrical pores. The denominator ensures volume conservation during PSD updating.

The total porosity after N cycles is obtained by integrating the updated PSD:

nN=∫φN(D)dD(20)

Using the initial PSD (e.g., Fig. 3), the model recursively computes the PSD for successive cycles (Fig. 6), capturing the migration of pore volume toward larger sizes and the corresponding increase in total porosity.

Figure 6: The aperture distribution curves after different number of F-T cycles

Fig. 6 illustrates the PSD curves predicted by the microscopic model for different numbers of F-T cycles. The simulation results align well with experimental trends: a substantial increase in the volume fraction of larger pores [2]. This shift visually demonstrates the cumulative effect of the non-recoverable frost deformation (1−Re(Di)), which is amplified to varying degrees in pores of different sizes. Thus, the microscopic model provides a mechanistic explanation for the macroscopic growth of porosity.

3.2 Effective Permeability Coefficient and Effective Pore Pressure

Within the multiphysics-coupled F-T damage model of concrete, the dynamic evolution of permeability is synergistically regulated by pore structure reconfiguration and phase transition processes. Based on pore-wall adsorption effects, the intrinsic permeability ki for individual pores with diameter Di can be formulated as:

ki=η(Di)kmaxDiDmax(21)

where the pore diameter correction coefficient η(Di) is assigned according to pore typology: η∈[0.3,0.5] for gel pores, η∈[0.5,0.8] for capillary pores, and η=1 for non-capillary pores [43]. kmax corresponds to the permeability at the maximum pore diameter Dmax.

The macroscopic permeability k0 is determined by integrating contributions from all pore sizes through weighted summation of the PSD function φ(Di):

k0=∫DminDmaxk0iφ0(Di)n0dDi=∫DminDmaxη(Di)kmaxDiDmaxφ0(Di)n0dDi(22)

in the equation, n0 represents the macroscopic initial porosity of the material. After Nth F-T cycle, the permeability coefficient kN decays due to cumulative porosity expansion (Eq. (20)). Consequently, kmax can be inversely derived as:

kmax=kN−1(∫DminDmaxη(Di)DiDmaxφN−1(Di)nN−1dDi)−1(23)

The effective permeability coefficient of unfrozen water kef is governed by the coupled effects of PSD function φ(Di) and phase transition thresholds Dcr. Given that the PSD function exhibits approximate Gaussian characteristics in logarithmic space, direct integration over linear scales induces numerical instability due to spanning multiple orders of magnitude. To address this, the integration kernel is transformed into logarithmic coordinates and expanded via Gauss-Legendre quadrature:

kef=∫DminDcrη(Di)kmaxDiDmaxφN−1(Di)nNdDi≈kmaxDmaxnN∑m=1M−1hm∑j=1Nωj(m)η(Dj(m))Dj(m)φN−1(Dj(m))(24)

Specifically, the interval [Dmin,Dcr] is discretized into M subintervals, with endpoints denoted as Dm and Dm+1. The half-length of each subinterval is defined as hm=(Dm+1−Dm)/2, and the midpoint as cm=(Dm+1+Dm)/2. Each subinterval is mapped to the physical pore-size scale via a linear transformation of the standard Gauss-Legendre nodes gj∈[−1,1]. For the mth subinterval, the jth Gauss node corresponds to a pore diameter of Dj(m)=cm+hmgi, where the nodal positions gj and corresponding weights ωj(m) are determined by the orthogonality properties of Legendre polynomials.

As evidenced in Fig. 7, under constant initial porosity conditions, the permeability coefficient exhibits a characteristic exponential decay during individual F-T cycles, due to increasing ice content as temperature reduction. However, with increasing cycle numbers, the initial permeability at the start of each cycle displays a distinct stepwise escalation. Notably, after multiple cycles, the temperature sensitivity of the permeability coefficient intensifies, manifesting as accelerated decay rates within identical temperature reduction intervals as cycling progresses. This phenomenon indicates that repeated F-T processes progressively damage the pore structure, increasing the proportion of large-diameter pores and rendering the overall pore network more fragmented with enhanced connectivity, thereby elevating the baseline permeability at the onset of subsequent cycles.

Figure 7: The permeability evolution with temperature at different FTC

Regarding the evolution of pore pressure during freezing, it is essential to distinguish between the total pore pressure and the effective pore pressure. The total pore pressure refers to the pressure within the pore fluid mixture. In contrast, the effective pore pressure is the net pressure acting on the solid skeleton, which is critical for calculating the effective stress that governs mechanical deformation.

Based on the Young-Laplace thermodynamic equilibrium equation [8], the pressure difference across the ice-water interface in fully frozen pores is expressed as:

Δp=pc−pL=ΥcL(1D1+1D2)(25)

where Δp denotes the pressure differential across the interface, pc denotes the ice crystallization pressure, pL denotes the permeation pressure. ΥcL represents the solid-liquid interfacial energy of ice-water interface. For spherical pores (D1=D2=Dcr), incorporating a thickness correction for the unfrozen water film, the ice-induced pressure exerted on pore walls becomes:

pc=pL+ΥcL(2Dcr−1Di−δ)(26)

During pore freezing, the effective pore pressure pef comprises two components: hydrostatic pressure from unfrozen water migration and crystallization pressure from ice formation. Consequently, the temperature-dependent evolution of pef is characterized as [50]:

pef(θ)=1nN∫DminDcrpL(θ)φ(Di)dDi+1n0∫DcrDmaxpδφ(Di)dDi=pL(θ)+ΥcLnN∫DcrDmax(2Dcr−1Di−δ)φ(Di)dDi(27)

where θ=Tf−T0(°C) denotes the degree of supercooling. And the second right item is defined as X(θ). Thus, the effective pore pressure is further simplified to:

pef(θ)=pL(θ)+X(θ)(28)

The evolution of the pressure parameter X(θ) is governed by the synergistic interplay between the phase-change threshold Dcr and the PSD curve φ(Di). Numerical implementation of Equation X(θ) via Gaussian quadrature quantitatively resolves the temperature-dependent characteristics of X(θ) under varying F-T cycles, as illustrated in Fig. 8. results demonstrate that progressive activation of phase transitions during cooling induces a monotonic increase in X(θ). Concurrently, cumulative pore-structure expansion from repeated cycling amplifies the thermal sensitivity of X(θ), yielding elevated pressure generation under identical cooling magnitudes.

Figure 8: Evolution of pressure parameter X(θ) during freezing

3.3 Governing Equation for THM System

3.3.1 Mechanical Equilibrium

The constitutive behavior of porous materials under F-T conditions can be derived from poroelasticity theory. The relationship between effective stress σe and total Cauchy stress σ is expressed as:

σ=σe−αBpefI(29)

where I represents the 2nd-rank identity tensor, αB is Biot coefficient [12] defined as αB=1−KKs, with K=E3(1−2v) and Ks denoting the bulk moduli of the porous system and solid skeleton, respectively.

Incorporating thermal expansion effects through the Duhamel-Neumann formulation [51], the constitutive relationship becomes:

σe=E0:εef=E0:(ε−εT)(30)

here ε,εef, and εT represent the total strain tensor, effective strain tensor, and thermal strain tensor, respectively. The thermal strain is given by εT=αs(T−Tref)I, where αs is the linear thermal expansion coefficient, T is the current temperature, and Tref denotes the reference temperature. The stiffness tensor E0 is decomposed into active E0+(λ,v) and inactive E0−(λ,v) components based on phase-field damage theory.

The mechanical equilibrium equation under static loading conditions is:

∇⋅σ+b=0(31)

where b represents body forces, which are negligible compared to hydraulic and crystallization pressures during freezing [43].

Substituting Eqs. (6), (29) and (30) into the equilibrium equation yields:

∇⋅{[g(ϕ)E0+(λ,v)+E0−(λ,v)]:(ε−αs(T−Tref)I)−αBpefI}=0(32)

with corresponding boundary conditions:

u(x,t)=u¯(x,t)  x∈∂Ωuσ(x,t)⋅nf¯(x)=0x∈∂Ωf¯(33)

This formulation comprehensively characterizes the coupled interactions among seepage, thermal, damage phase, and stress fields [23].

3.3.2 Mass Balance for Pore Water

Based on Powers’s static water pressure hypothesis [6], the mass variation of water within closed pores in a saturated porous medium results from the combined effects of water consumption due to ice formation and the mass flux of unfrozen water migrating across the boundary. The phase transition and water flux obey liquid mass conservation, expressed as:

d(nNscρc)dt+∇⋅(nNρLsLvL)=−d(nNsLρL)dt(34)

where ρL and ρc denote water and ice densities, and vL is the velocity of unfrozen water. The right-hand item of the equation describes water mass change, while the left-hand side accounts for ice mass change and water flux item.

Following Darcy’s law, pore water velocity relates to hydraulic head [52]:

nNsLvL=−kefμL∇pL(35)

where μL denotes the dynamic viscosity coefficient of water.

The temperature-dependent densities of ice and water during phase transition are explicitly incorporated [22,50]:

dρJdt=ρJ0(p˙JKJ−3αJT˙)J∈{c,L}(36)

Poromechanics theory [51,53] establishes quantitative relationships between pore structure evolution and F-T-induced damage through strain-porosity coupling equations, while incorporating cumulative porosity damage over successive F-T cycles. The governing equation is formulated as:

dnNdt=αB−nNKsp˙ef−3(αB−nN)αsT˙+(αB−nN)∂εvol∂t(37)

Combining these relationships yields the seepage field governing equation:

Ap˙L+BT˙+C∂εvol∂t−ρLkefμLΔpL=0(38)

where Δ is the Laplacian operator, and the coefficients A, B and C are:

A=nN(sLρL/KL+scρc/Kc)+(scρc+(1−sc)ρL)(αB−nN)/KsB=−3(scρc+(1−sc)ρL)(αB−n0)αs−3nN(sLρLαL+scρcαc)+[(scρc+(1−sc)ρL)(αB−nN)/Ks]∂X(θ)/∂TC=(scρc+(1−sc)ρL)(αB−nN)(39)

Material damage significantly alters hydraulic properties. The computational domain is divided into a reservoir zone Ωr(t) (intact region), fractured zone Ωf(t) (fully damaged region), and transition zone Ωt(t). Hydraulic parameters within the transition zone are determined through linear interpolation using indicator functions χr and χf, satisfying [54]:

{χr(⋅,ϕ)=c2−ϕc2−c1c1<ϕ<c2χf(⋅,ϕ)=ϕ−c1c2−c1c1<ϕ<c2(40)

where c1,c2 are two phase field thresholds. In this paper, c1,c2 are set 0.1 and 1.0 separately.

The hydraulic characteristics are then expressed as:

{k=kefχr+kfχfαB=αBrχr+αBfχfn=nNχr+nfχf(41)

where the subscripts r and f represent parameter values in reservoir and fractured zones, respectively. This study defines the initial permeability coefficient k0=8.32×10−21, the permeability coefficient in the fractured zone as kf=1×10−15. The Biot coefficients are defined as αBr=0.65,αBf=1, respectively.

The boundary conditions are specified as:

pL(t)=p¯L on ∂Ωp(nNnLρLvL)⋅n=qL on ∂Ωqp(42)

3.3.3 Energy Conservation for Heat Transfer

The temperature field in saturated porous materials is governed by the energy conservation equation:

(ρC)ef∂T∂t−∇⋅(λef∇T)−∂Q∂t=0(43)

where (ρC)ef represents the effective volumetric heat capacity, accounting for difference in heat capacity and thermal conductivity between ice, water, and matrix. This parameter is defined as [53]:

(ρC)ef=(1−nN)ρsCs+nNnLρLCL+nNncρcCc(44)

here ρs denotes matrix density, while Cs,CL,Cc represent the specific heat capacities of matrix, water and ice, respectively.

The effective thermal conductivity λef for saturated porous materials is expressed as [55]:

λef=λs+nN(1−nN)/3λs+λs/(sLλL+scλc−λs)(45)

where λs,λL,λc correspond to the thermal conductivities of the matrix, water, and ice.

The released heat rate ∂Q∂t is given by: ∂Q∂t=L∂mc∂t=L∂(nNscρc)∂t, where mc represents ice mass and L=333.5 kJ/kg denotes the latent heat of fusion [23,53].

The temperature-phase field coupling is established through the temperature-dependent thermal conductivity [56]:

λs(ϕ)=(1−ϕ)2λs0(46)

where λs0 is the initial matrix thermal conductivity.

The boundary conditions are specified as:

T(t)=T¯ on ∂ΩT(−λef∇T)⋅n=qT+h(T−Ta) on ∂ΩqT(47)

where qT represents the heat flux normal to the surface, and Ta denotes the environment air temperature [53].

These equations capture the coupled THM processes in the concrete during F-T cycles, allowing for the simulation of damage evolution and material behavior under such conditions.

3.4 Cross-Scale Coupling Modelling Scheme

This model establishes a comprehensive multiphysics framework to capture the key coupling pathways governing F-T damage evolution in saturated porous media. The phase-field variable interacts with the mechanical properties by degrading the elastic energy and simultaneously modulates the spatial distribution of permeability and thermal conductivity, thereby linking damage evolution with both seepage and thermal responses. The temperature field, influenced by ice-water phase transition and latent heat exchange, alters the pore-scale ice content, reshaping the pressure distribution and inducing thermal strains that feed back into the mechanical behavior. A bidirectional coupling exists between seepage and mechanical fields: pore pressure generated by freezing acts as an internal stress source, while deformation-induced volumetric changes affect pore connectivity and fluid transport. Stress redistribution drives the progression of damage through the phase-field, compounded by fatigue accumulation under cyclic loading. Repeated F-T cycles further modify the pore structure and permeability in a cumulative manner, forming a feedback loop that accelerates the degradation of transport properties and structural integrity, as illustrated in Fig. 9.

Figure 9: Cross-scale coupling modelling scheme

4  Numerical Discretization and Implementation

This section establishes a computational framework for modeling multiphysics-coupled fracture evolution in concrete under F-T cycles based on phase-field theory. The numerical implementation employs a hybrid space-time discretization strategy: the spatial domain is discretized using the finite element method, while temporal discretization is achieved via the backward differentiation formula. To effectively address numerical instabilities arising from strongly coupled multiphysics interactions, we propose a staggered iterative solution scheme coupled with an adaptive time-stepping control mechanism to balance computational efficiency and convergence robustness.

4.1 Weak Form and Discretization

The weak forms of the governing equations for F-T damage are expressed as:

∫Ω[2l0(1−kζ)HGc+f(α¯(t))]ϕδϕdΩ−l02∫Ω∇f(α¯(t))⋅∇ϕ⋅δϕdΩ−l02∫Ωf(α¯(t))⋅∇ϕ⋅∇δϕdΩ−∫Ω2l0(1−kζ)HGcδϕdΩ=0(48a)

−∫ΩDε:δεdΩ+∫ΩDαs(T−Tref)I:δεdΩ+∫ΩαBpefI:δεdΩ+∫∂Ωtf¯⋅δudS+∫ΩF:δudΩ=0(48b)

∫ΩAp˙LδpdΩ+∫ΩBT˙δpdΩ+∫ΩC∂εvol∂tδpdΩ+∫ΩρLkefμL∇pL∇δpdΩ−∫∂ΩtqLδpdS=0(48c)

∫Ω(ρC)ef∂T∂tδTdΩ+∫Ωλef∇T∇δTdΩ−∫ΩL∂mc∂tδTdΩ−∫∂Ωt(qT+h(T−T∞))δTdS=0(48d)

The spatial discretization scheme is constructed using the Bubnov-Galerkin weighted residual method [41], where the continuous field variables and admissible variations are formulated as:

{ϕ(x)=Nϕϕ~ δϕ(x)=Nϕδϕ~u(x)=Nuu~ δu(x)=Nuδu~ pL(x)=Npp~L δpL(x)=Npδp~LT(x)=NTT~ δT(x)=NTδT~(49)

the continuous field variables ϕ,u,pL,T and their variational counterparts δϕ~, δu~, δp~L, δT~ are spatially projected through linear combinations of shape function matrices Nϕ,Nu,Np,NT and their associated nodal degree-of-freedom vectors. This discretization rigorously preserves the Bubnov-Galerkin symmetry condition, wherein trial and test functions reside in identical Hilbert spaces.

Nu=[N100⋯Nn000N10⋯0Nn000N1⋯00Nn],Nϕ=Np=NT=[N1N2⋯Nn](50)

with Ni representing the shape function of the ith node in element.

Following analogous methodology, the gradient terms and corresponding variational formulations for each physical field are derived as:

{∇ϕ(x)=Bϕϕ~ ∇δϕ(x)=Bϕδϕ~ε(x)=Buu~ δε(x)=Buδu~∇T(x)=BTT~ ∇δT(x)=BTδT~∇pL(x)=Bpp~L ∇δpL(x)=Bpδp~L(51)

the differential operator matrices Bϕ,Bu,Bp,BT are the derivative of the shape functions, which can be derived as follow:

Bu=[N1,x00⋯Nn,x000N1,y0⋯0Nn,y000N1,z⋯00Nn,z0N1,zN1,y⋯0Nn,zNn,yN1,z0N1,x⋯Nn,z0Nn,xN1,yN1,x0⋯Nn,yNn,x0],Bϕ=Bp=BT=[N1,xN2,x⋯Nn,xN1,yN2,y⋯Nn,yN1,zN2,z⋯Nn,z](52)

To ensure the well-posedness of the discrete system, the selection of weight functions preserves mathematical arbitrariness, leading to the matrix formulation of the discretized governing:

[KϕKuSp+MpST+MT]{ϕ~u~p~LT~}−{FϕextFuextFpLextFText}={R~ϕR~uR~pLR~T}(53)

where R~ϕ, R~u, R~p, R~T represent the residual vector for equilibrium of fracture phase evolution flow, stress, seepage and thermal flow, respectively; Ku,Kϕ denote the stiffness matrices for the stress field and phase-field, respectively; Sp,ST correspond to the compression coefficient matrices and the phase transition heat capacity matrices, Mp,MT are the permeability matrix and the thermal conductivity matrix, Fϕext,Fuext,Fpext,FText signify the external force vectors for coupled fields.

4.2 Time Integration and Nonlinear Iterative Solution

The F-T cycling process inherently constitutes a complex physical system involving coupled THM phase transition interactions, exhibiting pronounced nonlinear characteristics and stiffness evolution over time. The ice-water migration and pore structure reconfiguration during phase transitions induce abrupt temporal gradient variations in field variables, imposing stringent requirements on the numerical stability and computational accuracy of temporal integration algorithms.

This study employs the backward differentiation formula for temporal discretization, which maintains numerical stability while achieving second-order temporal accuracy, making it particularly suitable for strongly nonlinear coupled systems. The discretized form of the temporal derivative term is expressed as:

∂x~∂t|k+1≈3x~k+1−4x~k+x~k−12Δt(54)

the variables x~(k−1),x~k and x~(k+1) represent the numerical solutions x at time steps k−1,k and k+1, respectively, with Δt denoting the time step size.

The updating algorithm for fatigue damage evolution parameters adopts an incremental accumulation scheme, where historical damage variables are progressively updated through cyclic integration of localized energy dissipation metrics.

α¯k+1=α¯k+∫tk+1tk+1α¯˙dτ=α¯k+Δα¯(55)

with Δα¯ can be approximated as

Δα¯=|αk+1−αk|H(αk+1−αkΔt)(56)

where H(⋅) denotes the Heaviside step function, enforcing the irreversibility of damage accumulation.

4.3 Solution Algorithm

This study develops a multiphysics-coupled numerical model on the COMSOL Multiphasics platform [41], employing a modular equation framework to simulate the complete F-T damage evolution process. The model incorporates six sets of core governing equations, as detailed in Table 2.

The computational framework employs a hybrid discretization strategy:

(1)   Spatial discretization: Phase-field variable ϕ, temperature T, seepage pressure pL, and displacement u are discretized using second-order Lagrange elements via COMSOL’s Partial Differential Equationmodule.

(2)   Temporal evolution: History-dependent variables H and fatigue damage α¯ are resolved through Ordinary Differential Equation module of COMSOL.

(3)   Nonlinear solver configuration: The system exhibits strong nonlinearity due to the stiffness degradation caused by ∂σe/∂ε and phase field. To ensure computational accuracy and stability, we employ a Newton-Raphson-like iterative scheme (NRI) [57] with a relative tolerance of 10−4 and absolute tolerance of 10−6 for all field variables. The staggered iterative scheme employs sequential solution of thermal, hydraulic, mechanical, and phase-field subproblems within each time step, with convergence achieved when the maximum relative change between successive iterations falls below 0.1%. Anderson acceleration enhances convergence for ill-conditioned Jacobians.

To implement the computational procedure effectively, a step-by-step algorithm is established as summarized in Algorithm 1, which outlines the sequential solving strategy for the F-T cycles.

5  Validation and Discussion

This section validates the proposed multiscale framework using experimental data from Wang et al. [58], which involve uniaxial compression tests on concrete subjected to rapid F-T cycles [59]. Two-dimensional symmetric model is constructed to replicate the test specimen, incorporating a mesostructured composed of aggregates, ITZ and mortar matrix. The accuracy and reliability of the proposed model will be verified by comparing the simulated temperature field and mechanical response curves with experimental measurements.

5.1 Geometrical Modeling and Parametric Framework

Utilizing the mirror symmetry characteristics of the specimen in orthogonal planes, this study establishes a quarter-symmetric model (50×50 mm2) to reduce computational costs and optimize efficiency. Symmetric constraints are enforced on orthogonal planes. Aggregates are generated as randomly shaped polygons via the Voronoi algorithm, maintaining statistically uniform spatial distribution in boundary regions without anomalous clustering along symmetry planes. Consistent with typical concrete microstructure characteristics, the thickness of ITZ between aggregates and the matrix is set within 0.03∼0.06 mm, which aligns with established experimental obsevations of interfacial transition zones in cementitious composites,as illustrated in Fig. 10.

Figure 10: The concrete sample with randomly distributed aggregates

To characterize the pore heterogeneity in the matrix material, a multi-scale self-similar porosity random field is generated using the diamond-square algorithm [60], as shown in Fig. 11.

Figure 11: Randomly distributed porosity

This modeling method captures porosity heterogeneity from the specimen scale down to the microscale, offering a high-resolution, statistically self-similar basis for simulating the spatially nonuniform evolution of F-T damage.

In order to enhance clarity, the material parameters are organized into three separate tables according to their physical roles in the model.

The phase-field fracture parameters associated with mortar, aggrates, and the ITZ are summarized in Table 3. Mortar and ITZ are modeled as temperature-dependent poroelastic media, while aggregates-due to their significantly higher stiffness and negligible permeability-are treated as non-damageable inclusions. The critical energy release is obtained through sensitivity analyses and previous investigations [61,62], and the phase-field length scale is set with reference to the characteristic material dimensions.

The THM parameters governing the multiphase F-T behavior of the three constituents are provided in Table 4. Mortar and ITZ properties are derived from experimental datasets and established literature, with interfacial properties further determined by Nagai et al.’s regression formula [63] based on water-to-cement ratios. Aggregate properties are simplified to thermal conduction only, given their impermeability. Fatigue-related parameter αT is calibrated against experimental observations to capture progressive deterioration under cyclic freezing and thawing, defined as: αT=5.626×10−1(N/mm).

Table 5 summarizes the material constants associated with the ice-water phase transition for F-T analysis [22,55,64]. The dynamic viscosity of water is defined as a temperature-dependent function, defined as [64]: μL=2.88×10−5exp⁡(509.53/(123.15+T)).

5.2 Boundary Conditions Configuration

Following ASTM C666-A standard procedure, specimens underwent 28-day water curing prior to F-T cycling in sealed rubber containers. Programmable thermal controllers imposed boundary conditions with cyclic protocols: each 6-h cycle included linear cooling from 20°C to −25°C at 0.75°C/min (60 min), 120-min isothermal holding, and reheating to initial temperature at identical rates (Fig. 12a). High-precision thermocouples and strain gauges were installed at the specimen surface (Point A) and geometric center (Point B), respectively, to monitor real-time temperature and strain evolution.

Figure 12: FTC boundary setting: (a) Cyclic temperature; (b) Boundary conditions

The boundary conditions in the numerical model were rigorously configured to replicate the experimental conditions (Fig. 12b). The surface temperature loading was applied using the measured temperature profile, while the temperature variation at the core was validated by comparing simulation results against the recorded data at monitoring Point B.

Mechanical boundary conditions assumed zero initial stress and strain states. Hydraulically, initial pore pressure was set to pL0=0, with zero boundary pressure maintained throughout cycling, preventing external fluid exchange.

Following 75 F-T cycles, specimens were subjected to 48-h isothermal curing in a temperature-controlled water tank (20°C ± 1°C) to dissipate residual thermal stresses, prior to uniaxial compression testing. The compressive damage simulation employed the following boundary configurations: the bottom surface of the specimen was fully constrained in all degrees of freedom, while the top surface was subjected to displacement-controlled loading with a step size of 0.1 mm, as illustrated in Fig. 13. The material parameters in the numerical model were assigned based on the experimentally measured elastic modulus and Poisson’s ratio of the F-T damaged specimen to accurately capture the cumulative damage effects.

Figure 13: Simulation model for compression

5.3 F-T Damage and Static Strength

Fig. 14 illustrates a comparative analysis of the temperature evolution at the specimen’s geometric center (Point B) between experimental measurements and numerical simulations during a single F-T cycle. The simulated results exhibit consistent trends with experimental data, accurately capturing the phase transition plateau near 0°C caused by latent heat release during water-ice transformation. The plateau duration is governed by the material’s thermal conductivity, pore structure, and saturation conditions.

Figure 14: Temperature variation curve of material center point (B)

To introduce varying degrees of F-T damage prior to analyzing axial compression responses, unconstrained expansion during 75 F-T cycles (FTCs) was simulated, corresponding to experimental procedures by Wang et al. [58]. Fig. 15 presents simulated results of free FTC expansion, serving as initial conditions for subsequent axial compression analysis. The residual strain after each FTC constitutes the most critical parameter, which can be linked to other material properties [65]. Consequently, once simulated residual strains match experimental measurements, FTC-induced damage observed in tests can be replicated by the model. Fig. 15a compares simulated and experimentally measured FTC strain histories, demonstrating that both datasets exhibit consistent net expansive strain accumulation, with simulated trends aligning closely with experimental observations. Quantitative analysis reveals a post-75-cycles residual strain of 1104 µε (experimental) vs. 1057 µε (simulated), yielding a 4.2% relative error-demonstrating the model’s engineering reliability in predicting long-term damage accumulation.

Figure 15: Strain history comparison: (a) Simulated cyclic strain and experimental measurements in 75 FTC; (b) First cycle strain

The dynamic strain response during the initial F-T cycle displays characteristic biphasic behavior in numerical results as illustrated in Fig. 15b. This behavior arises from competing thermal and phase-change effects: during initial cooling stages, rapid surface thermal contraction generates compressive strains near the boundaries. However, as temperatures continue to decrease below the critical freezing point, pore water begins to freeze and expand, generating significant crystallization pressure that reverse the strain direction and produce net expansion. The transition occurs when ice-induced expansion forces overcome thermal contraction effects [58].

The framework effectively resolves transient interactions between thermal gradients and phase-change-induced stresses, capturing localized damage initiation at material heterogeneity zones. Crack propagation pathways demonstrate preferential alignment with interfacial transition regions, consistent with microscopic observations of frost-damaged specimens. This validation confirms the model’s capacity to simulate coupled thermo-mechanical degradation processes while maintaining physical consistency across scales.

Fig. 16 compares the axial compressive failure modes between pristine specimen (Fig. 16a) and one subjected to 75 F-T cycles (Fig. 16b). Numerical simulations reveal that the failure of undamaged specimens is dominated by shear bands, exhibiting typical X-shaped intersecting crack propagation paths. In contrast, F-T damaged specimens demonstrate distributed network-like crack growth along interfacial transition zones (ITZ), with cracks preferentially developing at weakened aggregate-mortar interfaces, the simulation results are consistent with experimental observations [58]. This phenomenon originates from the pre-damaged effect induced by F-T cycles in ITZ regions, repeated phase transitions generate microcracks at interfaces, significantly degrading their mechanical properties. Consequently, subsequent compressive loading concentrates damage in pre-damaged areas, creating a synergistic amplification mechanism between F-T and mechanical degradation. Both experimental and numerical results uncover the reconfiguration of material failure modes by frost damage: F-T pre-damage not only reduces overall strength but also alters crack propagation pathways, thereby modifying failure mechanisms. The numerical model successfully reproduces the transition from shear-dominated failure mechanisms. This phenomenon reveals the reconfiguration mechanism by which F-T damage alters material failure modes, providing a theoretical foundation for predicting multi-stage coupled degradation processes.

Figure 16: Damage patterns under axial compression: (a) FTC_0; (b) FTC_75

Fig. 17 compares the experimental (dashed line) and simulated (solid line) stress-strain curves of concrete specimens under axial compression, normalized against the room-temperature peak strain (εpeak≈0.002) [16,23]. Extract the compressive strengths and simulation results across varying F-T cycles from the Fig. 17. After 75cycles, experimental and simulated peak stresses decrease to 67% and 65% of initial values, respectively, validating the significant strength degradation caused by F-T action. The post-peak softening stages of cycled specimens exhibit pronounced prolongation, reflecting internal crack development and stress redistribution during F-T processes [58]. Notably, potential overestimation of peak strain in experimental measurements due to testing apparatus constraints has been systematically accounted for in result analysis.

Figure 17: Compressive stress-strain curves

The simulation accurately captures two key phenomena: (1) strength reduction and ductility enhancement, and (2) stiffness recovery at higher stress levels. Analysis of the phase-field variable evolution during compression reveals that this stiffness recovery primarily results from microcrack closure under increasing compressive loads. As stress levels rise, existing microcracks gradually close, restoring some load transfer capacity through improved interfacial contact. Additionally, stress redistribution occurs from heavily damaged regions to relatively intact areas, further contributing to the apparaent stiffness recovery. This mechanism is particularly pronounced in F-T damaged specimens dur to the distributed nature of frost-induced microcracking. This partial stiffness recovery, emerging after initial degradation, provides new insights into concrete’s performance restoration mechanisms following F-T damage.

5.4 Discussion

5.4.1 The Evolution of F-T Damage

The multiphysics-coupled model enables controlled characterization of frost damage progression, as illustrated by damage contour maps after 25, 50, 75, and 100 F-T cycles (Fig. 18a–d). Numerical investigations reveal that initial damage nucleation preferentially occurs in regions of significant microstructural heterogeneity, driven by thermally induced pore water migration and nonuniform stress fields:

Figure 18: Comparison of morphology and simulation result of concrete F-T damage: (a) FTC = 25; (b) FTC = 50; (c) FTC = 75; (d) FTC = 100; (e) Experimental observatiaon

(1)   Surface damage initiation: Thermal-hydraulic-mechanical interactions concentrate volumetric strain from ice-water phse change in the weakest boundary zones, leading to nucleation of microcracks at the surface;

(2)   ITZ-dominated propagation: Stress concentration and crack advancement preferentially follow ITZ due to dual weakening effects-reduced elastic modulus and compromised compactness.

As cycling continues, damage spreads progressively. Within the ITZ, capillary migration of unfrozen water and tensile stresses from ice crystal growth synergistically drive cracks from the aggregate-matrix interface toward the specimen center. These microcracks gradually interconnect, producing crack morphologies that closely match experimental observations (Fig. 18e). Concurrently, surface initiated cracks coalesce, ultimately resulting in spalling and surface delamination.

After 100 F-T cycles, the simulated damage pattern-dominated by crack networks along the ITZ-exhibits strong spatial correlation with microscopic imagery, confirming the model’s capability to capture the underlying mesoscopic damage mechanisms under coupled thermal, hydraulic, and mechanical loading.

5.4.2 F-T Hysteresis Phenomenon

Fig. 19 illustrates the spatiotemporal evolution of ice saturation at boundary Point A and central Point B during F-T cycles, along with their thermodynamic path dependence, characterized as follows:

Figure 19: Ice saturation evolution: (a) Temporal evolution of ice saturation; (b) Temperature-dependent evolution of ice saturation

(1)   Temporal evolution (Fig. 19a): The peak ice saturation at Point B lags significantly behind Point A, with a 4% reduction in peak magnitude. This delay reflects spatiotemporal heterogeneity in heat conduction under temperature gradients, where slower thermal equilibration in the interior governs phase-change dynamics.

(2)   Temperature-dependent evolution (Fig. 19b): Ice saturation-temperature curves form closed hysteresis loops, with freezing-phase supercooling exceeding thawing-phase values. This asymmetry aligns with phase-transition hysteresis mechanisms documented in [46], confirming the thermodynamic irreversibility of ice-water transformations.

Faster heat exchange at boundary regions creates divergent phase-change kinetics compared to the interior. Such spatiotemporal asynchrony generates non-equilibrium stress fields that redirect damage evolution paths, explaining the gradient-driven accumulation of F-T damage.

5.4.3 The Influence of Water-Cement Ratio on F-T Damage

The water-cement (w/c) ratio critically governs the microstructure andmacroscopic mechanical properties of concrete, directly influencing frost resistance [63,66]. Adopting Kosaka et al.’s model [67] that quantifies w/c-dependent elastic modulus, specimens with w/c = 0.4, 0.5 and 0.6 were analyzed to assess damage sensitivity.

The relative elastic modulus loss rate (LE=ΔE/E0×100%) [68] quantifies frost damage progression. As shown in Fig. 20, LE exhibits nonlinear degradation:

Figure 20: Relative loss rate of Elastic modulus: (a) Bar chart of elastic modulus; (b) Line chart of elastic modulus loss rate

At 25 cycles, LE values for w/c = 0.4, 0.5, and 0.6 as 2.4%, 4.6%, and 7.8%, respectively, indicating that early-stage freezing and thawing induces only minor structural impairment, without substantially disrupting the internal load-transfer network. However, damage accumulates rapidly with additional cycling.

At 100 cycles, LE escalates to 25.9%, 30.3%, and 33.8% for the respective w/c ratios. This accelerated degradation is chiefly driven by tensile stresses from ice-water phase transitions within the pore space, which generate localized stress concentrations, nucleate microcracks, and promote their rapid propagation along interconnected pore channels. The resulting microcrack network precipitates a steep decline in the material’s load-bearing capacity.

Higher w/c ratios amplify damage due to dual mechanisms:

(1)   Elevated capillary porosity: Increased pore connectivity enhances frost heave pressure transmission;

(2)   Unhydrated particle accumulation: Residual unreacted particles reduce cementitious phase strength [63,66].

The synergy between porous network development and weak phase concentration accelerates crack propagation and structural deterioration under cyclic freezing-thawing.

5.4.4 Impact of Air Content on F-T Damage

As illustrated in Fig. 21, the introduction of entrained air bubbles into concrete significantly optimizes pore structure, enhancing frost resistance. Fig. 21a presents PSD curves for specimens with varying air contents. Increasing air content shifts the distribution from a broad scatter to a concentrated range of small-to-medium pores, indicating improved uniformity and controllability of the pore structure.

Figure 21: Average hydraulic pressure: (a) The PSD function [69]; (b) The effective permeability coefficient; (c) Hydraulic pressure

Based on the PSD functions in Fig. 21a, the temperature-dependent evolution of effective permeability coefficients for specimens with different air contents is calculated as illustrated in Fig. 21b. Results demonstrate that higher air contents correlate with denser pore networks and significantly reduced permeability.

Fig. 21c compares the temporal evolution of internal pore pressures during F-T cycles. Specimens with higher air contents consistently maintain lower average pressure levels. For instance, the maximum internal pressure in specimens with 5.6% air content is 12 MPa, representing a 52% reduction compare to non-air-entrained specimens (25 MPa). This highlights the buffer effect of air bubbles in mitigating ice-induced pressure peaks, thereby delaying crack initiation and propagation. The trend suggests that increased air content shortens capillary water migration paths during freezing, reducing hydrostatic pressure generation [6,12] and alleviating localized damage.

In summary, optimizing air content enables synergistic improvements in pore structure and permeability under multiphysics coupling, effectively suppressing peak hydraulic pressures and significantly enhancing concrete’s frost durability.

6  Conclusions

This study proposes a novel multiphysics-coupled computational framework for F-T damage, with core innovations extending existing phase-field variational fatigue theory to thermos-hydro-mechanical systems and establishing cross-scale linkages between dynamic porosity evolution and fatigue damage. Specifically, leveraging the variational nature of phase-field methods, the damage phase-field variable is bidirectionally coupled with thermal conduction, hydraulic transport, and mechanical constitutive equations. Critically, the framework abandons conventional static porosity assumptions and introduces a cycles-dependent dynamic porosity evolution equation to iteratively characterize progressive microstructural degradation. This equation treats post-cycles porosity as the initial condition for subsequent cycles, inherently capturing damage accumulation path dependency. By mapping dynamic porosity to permeability evolution, the model quantitatively resolves transport property deterioration during F-T cycling.

Numerical investigations reveal key mechanisms:

(1)   Phase-transition dynamics: Initial F-T stages induce contraction from rapid cooling, followed by ice-expansion-dominated volumetric growth, confirming the critical transition mechanism from contraction to expansion.

(2)   Mesoscale damage evolution: During the initial stage (<25 cycles), discrete micro-damage localized near aggregate boundaries and ITZ. As cycling reaches 50, surface damage coalesces into connected networks, while internal damage remains localized. Beyond 75 cycles, internal microcracks propagate along ITZ, forming networked crack trends with exacerbated surface spalling.

(3)   Mechanical degradation: Post-F-T compressive strength decreases by 45%, with higher water-cement ratios accelerating performance deterioration, highlighting porosity’s pivotal role in frost susceptibility.

(4)   Air entrainment efficacy: Optimized pore structure via air entrainment reduces permeability and suppresses hydraulic pressures, diminishing peak damage stress attenuation by 18%-22%, validating microstructural control as an effective frost-resistance strategy.

The model successfully reproduces the transition from isolated microcracks to interconnected fracture networks while quantitatively predicting macroporosity evolution and its correlation with mechanical degradation. This capability provides a robust tool for assessing durability thresholds and optimizing frost-resistant concrete designs in cold regions. This study advances the state of the art by developing a comprehensive cross-scale modeling framework that effectively integrates dynamic pore evolution with Multiphysics coupling during F-T damage processes. Unlike previous studies that often treated pore structure as static or employed empirical degradation models, our approach resolves inherent constraints in conventional approaches regarding the characterization of damage path dependency and time-dependent microstructural evolution through a physics-based porosity evolution equation.

The proposed framework provides a general paradigm for simulating freeze-thaw failure. But it should be pointed out that the F-T evolution model involving pore structure is not yet refined enough, and still relies on phenomenological construction, without clearly reflecting the mechanism of pore structure thawing recovery. Future work will focus on establishing a more physics-based evolution law of pore structure.

Future work will focus on three main directions:

(1)   Multi-stressor coupling: Extending the model to account for combined environmental streesors, such as F-T cycles coupled with chloride ingress or salt crystallization, which exacerbate damage through synergistic chemical and mechanical effects.

(2)   Model optimization and validation: Conducting systematic sensitivity analyses of key material parameters and experimental validation under diverse environment conditions to enhance model robustness and applicability.

(3)   Computational efficiency improvements: Developing reduced-order modeling techniques and parallel computing strategies to enable larger-scale structural simulations while maintaining computational accuracy.

This phase-field-based framework not only advances fundamental understanding of F-T damage mechanisms but also provides a versatile platform for durability-oriented concrete design in cold climates. Subsequent research will prioritize translating these computational insights into field-applicable guidelines for infrastructure resilience.

Acknowledgement: The financial support of the National Natural Science Foundation of China is gratefully acknowledged.

Funding Statement: This work has been supported by the National Natural Science Foundation of China (Grant Nos. 11932006 and 12172121).

Author Contributions: The authors confirm contribution to the paper as follows: study conception and design: Siwei Zhang, Xiaozhou Xia, Qing Zhang; data collection: Siwei Zhang; analysis and interpretation of results: Siwei Zhang, Meilin Zong, Xin Gu; draft manuscript preparation: Siwei Zhang, Xiaozhou Xia. All authors reviewed the results and approved the final version of the manuscript.

Availability of Data and Materials: Readers can access the data used in the study by contacting email glzhangsiwei@hhu.edu.cn.

Ethics Approval: This study did not involve any human participants or animal subjects. As the research exclusively employed numerical simulations and theoretical modeling, ethical approval was not required for this work.

Conflicts of Interest: The authors declare no conflicts of interest to report regarding the present study.

References

1. Guo J, Sun W, Xu Y, Lin W, Jing W. Damage mechanism and modeling of concrete in freeze-thaw cycles: a review. Buildings. 2022 Aug 28;12(9):1317. doi:10.3390/buildings12091317. [Google Scholar] [CrossRef]

2. Zhou R, Mao H, Yang S, Wang Z, He X. Research on mechanical properties and pore structure evolution process of steam cured high strength concrete under freeze thaw cycles. Sci Rep. 2024 Nov 11;14(1):27530. doi:10.1038/s41598-024-76631-3. [Google Scholar] [PubMed] [CrossRef]

3. Chen C, Zhang K, Ye L. Influence of freeze-thaw cycles and sustained load on the durability and bearing capacity of reinforced concrete columns. Materials. 2024 Dec 15;17(24):6129. doi:10.3390/ma17246129. [Google Scholar] [PubMed] [CrossRef]

4. Zhang K, Zhou J, Yin Z. Experimental study on mechanical properties and pore structure deterioration of concrete under freeze-thaw cycles. Materials. 2021 Nov 1;14(21):6568. doi:10.3390/ma14216568. [Google Scholar] [PubMed] [CrossRef]

5. Kristensen PK, Niordson CF, Martínez-Pañeda E. An assessment of phase field fracture: crack initiation and growth. Phil Trans R Soc A. 2021 Aug 9;379(2203):20210021. doi:10.1098/rsta.2021.0021. [Google Scholar] [PubMed] [CrossRef]

6. Powers TC. A working hypothesis for further studies of frost resistance of concrete. J Am Concr Inst. 1945;16(4):245–72. doi:10.14359/8684. [Google Scholar] [CrossRef]

7. Powers TC, Helmuth RA. Theory of volume changes in hardened portland cement paste during freezing [Internet]. Skokie, IL, USA: Portland Cement Association; 1953(Bulletin (Portland Cement Association. Research and Development Laboratories. Research Department)). doi:10.14359/20128. [Google Scholar] [CrossRef]

8. Scherer GW. Crystallization in pores. Cem Concr Res. 1999 Aug;29(8):1347–58. doi:10.1016/S0008-8846(99)00002-2. [Google Scholar] [CrossRef]

9. Valenza II JJ, Scherer GW. Mechanism for salt scaling. J Am Ceram Soc. 2006;89(4):1161–79. doi:10.1111/j.1551-2916.2006.00913.x. [Google Scholar] [CrossRef]

10. Setzer M, Auberg R. Basis of testing the freeze-thaw resistance—surface and internal deterioration. In: Frost resistance of concrete. Boca Raton, FL, USA: CRC Press/E&FN Spon; 1997. p. 157–73. [Google Scholar]