1   Introduction

Soil water redistribution controls vegetation growth and regional hydrological processes. Seasonal freezing affects a large portion of the global land surface. In these regions, freeze-thaw cycles alter pore structure and redistribute water and solutes (Lv et al., 2022). These changes modify both hydraulic conductivity and thermal conditions. Leuther and Schlüter (2021) showed that freeze-thaw cycles can break soil aggregates and push water and solutes toward freezing fronts. As temperature decreases, pore water begins to freeze. This process reduces effective porosity and limits flow pathways. At the same time, some water remains unfrozen. This unfrozen water migrates under temperature gradients and matric suction (Reinmann et al., 2019). Water therefore accumulates near the freezing front, especially in shallow soil layers (Figure 1). Coastal frozen soils show additional complexity. In these areas, interactions among sea ice, permafrost, and groundwater are common. These interactions affect erosion rates and carbon transport. The behavior differs from inland frozen soils. Understanding these processes is important for infrastructure stability and environmental assessment.

Figure  1  Schematics of soil freezing processes and gain in shallow total soil water content (Li et al., 2025)

Download: Full-Size Img

Simulating water flow and heat transfer during freeze-thaw processes remains difficult. Conventional methods, such as the finite difference method (FDM), finite element method (FEM), and finite volume method (FVM), treat soil as a homogeneous continuum. These approaches operate at scales much larger than individual pores. They are suitable for large-scale groundwater flow. However, they cannot resolve pore-scale features. These include menisci, local connectivity loss, and structural changes caused by phase transitions. This limitation becomes more significant under freezing conditions. Water flow, heat transfer, and phase change are strongly coupled and highly nonlinear (Li and Yin, 2025). Direct use of macroscopic parameters often leads to large uncertainty. Soil heterogeneity and scale effects reduce the reliability of these parameters.

The lattice Boltzmann method (LBM) provides a mesoscale modeling approach for these problems. It evolves particle distribution functions on discrete lattices. This structure allows direct representation of complex pore geometry. It also captures fluid-fluid and fluid-solid interfaces as they evolve over time. LBM can handle multicomponent interactions without relying entirely on homogenized constitutive relations (Chen and Doolen, 1998; Han and Cundall, 2013). These features make it suitable for studying seepage, phase change, and freeze-thaw processes in porous media. In recent years, LBM has been applied to micropore-scale flow, unsaturated seepage, and phase-change heat transfer. It has also been used to investigate coupled thermo-hydraulic processes. These studies show that LBM can reveal pore-scale mechanisms that are difficult to capture with continuum models.

Research on seepage and freeze-thaw processes has increased in recent years. However, studies based on the lattice Boltzmann method remain limited and scattered. A clear synthesis is still lacking. In particular, the link between pore-scale mechanisms and larger-scale hydrological or engineering behavior is not well established. This issue is more evident in freezing environments and coastal regions, where boundary conditions and phase behavior are more complex. This paper reviews recent advances in LBM modeling of porous media. The focus includes saturated and unsaturated seepage, ice-water interface behavior, and freeze-thaw cycles. Three aspects are considered. These are model applicability, scale bridging and parameter interpretation, and the limitations of current approaches under complex environmental conditions.

2   Physical basis of seepage and freeze-thaw processes

2.1   Characteristics of pore structure

The soil is a heterogeneous porous medium whose seepage and freeze-thaw processes are highly regulated through the pore structure. The major structural features involved are porosity, pore size distribution, connectivity, and tortuosity that altogether define how much water can be stored, what kind of flow paths exist, and the resistance to fluid and heat transport. These features will have significant effects on phase-change temperature, ice distribution, and hydraulics near the freezing front in cold soils (Chen et al., 2015; Li and Yin, 2025; Zheng et al., 2016).

In this set of parameters, porosity refers to the relative volume fraction of pores, however, pore-size distribution and connectivity dictate if the pores will be able to contribute efficiently to seepage. Tortuosity denotes the complexity of flow paths and has a strong association with permeability and efficiency of mass transport. With soil pores at various levels, including smaller intra-aggregate pores and bigger inter-aggregate passages, there are several different pore domains that take on distinct functions regarding the retention of moisture, drainage, and water redistribution during freeze-thaw processes. Particularly, large pores are much more probable to act as major seepage channels while smaller pores tend to hold water and support unfrozen water when temperatures fall below zero degrees Celsius.

The complete characterization of soil pore geometry in numerical simulations has been difficult. The central point of concern concerning seepage and freeze-thaw analysis is not to copy all geometric details but to focus on those structures which play the strongest role in transport and phase transformation. Some previous research has found out that changes in micropore properties may impact local phase-change behavior, seepage characteristics at the ice front, and movement of water within unfrozen and frozen areas (Engemann et al., 2004; Saruya et al., 2014; Watanabe and Flury, 2008). This means that pore structure must be considered as one of the determining factors of the seepage evolution and freeze-thaw interaction since it has an immediate effect on effective fluid flow spaces and transport channels that decide large-scale seepage behavior.

2.2   Pore-scale seepage characteristics

Permeability is an important property of porous media, and the permeability coefficient is influenced by many factors. Among them, porosity, particle size, pore type, pore tortuosity, and particle-size distribution are major controlling factors. The seepage behavior of porous media generally follows Darcy’s law. Under low Reynolds number conditions, the average velocity of a saturated viscous fluid flowing through a porous medium is proportional to the pressure gradient ∂p/∂x_i and inversely proportional to the fluid viscosity μ(Carman, 1997):

where $v_{i}$ is the outlet velocity in the $i$-direction, $k_{i}$ is the intrinsic permeability coefficient, $K_{i}$ is the hydraulic permeability, $\mu$ is the dynamic viscosity, $p$ is the applied inlet pressure, and $x_{i}$ is the distance in the $i$-direction.

There are many approaches for predicting the intrinsic permeability of soils. One of the earliest and most fundamental formulations is the Kozeny-Carman (KC) equation proposed by Carman (1997), which is expressed as:

where $c$ denotes the Kozeny constant. Subsequently, Koponen et al. (1997) modified this equation by introducing the hydraulic tortuosity $T_{r}$ to describe the tortuosity of streamlines in porous media, yielding:

These relations show that permeability is not controlled by pore volume alone but is also influenced by pore connectivity and the tortuosity of flow paths. In pore-scale seepage analysis, effective porosity plays a key role, since only pores that are connected to continuous flow pathways contribute to transport. As total porosity decreases, the number of inactive pores increases, and flow behavior becomes more sensitive to connectivity loss and local geometric constraints. Permeability should therefore be treated as a result of pore structure, including connectivity, tortuosity, and effective flow space, rather than being described by porosity alone. This structural control becomes more pronounced under freezing and thawing conditions, where phase change and water redistribution continuously modify pore connectivity and seepage pathways (Ahmadi et al., 2011; Kaviany, 1995).

2.3   Physical mechanisms of freeze-thaw processes

Soil freeze-thaw cycles occur under non-equilibrium conditions. Heat conduction, phase change, and moisture migration are strongly coupled, and the system does not reach a steady state. As temperature decreases, pore water freezes and releases latent heat. This heat forms a temperature plateau near the freezing front and slows further ice growth. During warming, ice absorbs latent heat and melts. The temperature response is therefore delayed relative to external heating, which leads to different thermal behavior during freezing and thawing (Tian et al., 2014). Pore size also affects these processes. Ice tends to form first in larger pores and can block main drainage pathways. Smaller pores retain unfrozen water films, which allow limited hydraulic connectivity even at subzero temperatures.

Unfrozen water plays a key role in freeze-thaw seepage. Capillary forces, matric suction gradients, and temperature gradients drive this water toward freezing fronts according to Watanabe and Mizoguchi (2002) and Shi et al. (2024). This process supports ice lens growth and affects where ice accumulates and how pore structure changes. During thawing, these effects are partly reversed but do not fully disappear. According to Klöffel et al. (2024), ice melting reconnects pore pathways and can lead to rapid increases in permeability, while capillary conductivity recovers differently depending on the previous freeze-thaw history. This difference reflects a hysteresis effect and separates freeze-thaw seepage from isothermal flow behavior.

Repeated freeze-thaw cycles restructure soil at the pore scale. Leuther and Schlüter (2021) and Chuvilin et al. (2022) reported that aggregates can break, macropore fractions can change, and connected pathways can reform. These structural changes affect heat conduction and water migration. They also alter the coupling between thermal and hydraulic processes. Such behavior cannot be described accurately by a single macroscopic parameter. It therefore provides a basis for using pore-scale numerical methods, such as the LBM, to study these processes.

3   LBM models relevant to seepage and freeze-thaw processes

The LBM has been widely used to simulate pore-scale seepage, multiphase flow, and phase-change heat transfer in porous media. In soil-related multiphase simulations, several key factors need to be considered, including density contrast, surface tension, wettability, and phase-change coupling. The commonly used multiphase LBM models include the pseudopotential model, free-energy model, color-gradient model, and phase-field model. Compared with the pseudopotential model, the other three approaches require more complex treatments for wettability and phase-change coupling and are generally less efficient computationally. Existing studies suggest that improved versions of the Shan-Chen pseudopotential model, combined with the multiple-relaxation-time (MRT) scheme, realistic equations of state, and dynamic wettability treatment, perform well in simulating soil processes. In practical applications, this model can capture water-gas and water-ice interfaces without explicit interface tracking. It also supports relatively large density ratios and stable viscosity contrasts under the MRT framework. Wettability can be adjusted through fluid-solid interaction parameters to match contact angles. In addition, coupling with enthalpy-based thermal models allows the description of latent heat release and ice lens growth during freeze-thaw processes, while maintaining reasonable computational efficiency for pore-scale simulations.

3.1   Pseudopotential model

The pseudopotential model is one of the most widely used multiphase flow models within the LBM framework. The advantage of the SC approach is its intrinsic simplicity and mesoscopic nature. Its core concept is to introduce an additional force term into the LBM collision operator to mimic intermolecular interactions, thereby inducing phase separation and generating interfacial tension. Shan and Chen (1993) first proposed this model, commonly referred to as the Shan-Chen (SC) model. In this method, a pseudopotential is defined based on a weighted sum of particle densities at neighboring lattice nodes, which gives rise to a non-ideal interaction force that drives the fluid to separate into distinct phases (Shen et al., 2025). In early implementations of the pseudopotential model, the pseudopotential interaction force was incorporated into the lattice Boltzmann equation by modifying the equilibrium distribution function,

where $f_{\alpha}$ denotes the distribution function of phase $\alpha, \tau$ is the relaxation time, $f_{\alpha}^{\text {eq }}=\omega_{\alpha} \rho\left[1+\frac{\boldsymbol{e}_{\alpha} \cdot \boldsymbol{u}}{c_{s}^{2}}+\frac{\left(\boldsymbol{e}_{\alpha} \cdot \boldsymbol{u}\right)^{2}}{2 c_{s}^{4}}-\frac{u^{2}}{2 c_{s}^{2}}\right]$ is the equilibrium function, and the equilibrium velocity is defined as $\boldsymbol{u}^{\mathrm{eq}}=\left(\sum\limits_{\alpha=0}^{Q} f_{\alpha} \boldsymbol{e}_{\alpha}+\tau \delta_{t} \boldsymbol{F}_{\alpha}\right) /\rho$, where $\boldsymbol{F}_{\alpha}$ represents the pseudopotential interaction force used to describe intermolecular interactions within the fluid, which is given by:

where is the pseudopotential at spatial position $\boldsymbol{r}, \psi(\boldsymbol{r}+ \boldsymbol{e}_{\mathrm{i}} \delta_{t}$) is the pseudopotential at the neighboring node in the $\boldsymbol{e}_{\alpha}$ direction, $G$ denotes the interaction strength, and $\omega_{i}$ is the weighting coefficient (Shan, 2006). However, in soil seepage simulations, not only fluid-fluid interactions are involved, but also solid phases, etc.; therefore, the $\boldsymbol{F}_{\alpha}$ in $\boldsymbol{u}^{\mathrm{eq}}$ becomes:

where $\boldsymbol{F}_{\alpha}^{\text {coh }}, \boldsymbol{F}_{\alpha}^{\text {ads }}, \boldsymbol{F}_{\alpha}^{\text {b }}, \boldsymbol{F}_{\alpha}^{\text {T }}$ denote the cohesive force between fluids, the adhesive force between fluid and solid, the body force, and the phase-change driving force, respectively. The calculation procedure for each term is as follows (Martys and Chen, 1996):

where $\psi_{\alpha}$ is the interaction potential, $s_{i}$ is the floating variable used to determine the fluid state, with $s_{i}\left(\boldsymbol{r}+\boldsymbol{e}_{i} \delta_{t}\right)=0$ indicating a fluid node and $s_{i}\left(\boldsymbol{r}+\boldsymbol{e}_{i} \delta_{t}\right)=1$ indicating a solid node. $G^{\text {coh }}$ is the cohesion coefficient that regulates the surface tension between different components in binary fluid mixtures. $G^{\text {ads }}$ is the adhesion coefficient required to determine the wettability of different components, as well as the wetting characteristics. $G^{\mathrm{T}}, \rho_{\alpha}, \nabla \varphi$ are the temperature potential coefficient, soil density, and total water potential, respectively.

3.2   The LBM model for phase-change heat conduction problems

In frozen soil, moisture migrates along the temperature gradient from high to low temperatures during freezing and thawing processes; therefore, accurate numerical models for simulating temperature distribution in frozen soil are of critical importance. In LB models, the heat conduction process can be implemented by introducing an additional thermal lattice. Similar to the previous subsection, the LB equation for thermal diffusion based on the temperature distribution function can be expressed as follows (Shan, 1997):

To address phase-change heat conduction problems, Jiaung et al. (2001) proposed an LB model with a source term incorporated into the equation, which can be written as:

where $L$ is the latent heat coefficient, $C$ is the specific heat capacity, and $\varphi_{l, \alpha}$ is the liquid volume fraction. The purpose of introducing this term is to address the unsteady heat conduction process involving liquid-solid phase change in frozen soil, rather than ordinary heat transfer problems. For the classical heat conduction equation, the form without phase change is:

However, in simulating soil freezing and thawing processes, the latent heat term satisfying energy conservation must also be incorporated:

where $\frac{\partial \varphi_{l}}{\partial t}$ is the phase-change rate. The latent heat term cannot be neglected; otherwise, it will lead to erroneous freezing positions, non-physical abrupt changes in the temperature field, and severe distortion of the freezing front propagation velocity.

For seepage and freeze-thaw problems in porous media, the main value of thermal LBM formulations lies not only in solving the temperature field itself, but also in enabling coupled treatment of heat transfer, liquid-ice phase change, and interface migration within the same framework. This capability is especially important when freezing modifies pore connectivity, hydraulic behavior, and local seepage pathways. In practical applications, thermal formulations are therefore often combined with multiphase seepage models and enthalpy-based phase-change treatments to support the simulation of coupled water-heat-phase change processes.

4   The modeling framework and research progress of LBM in seepage and freeze-thaw processes

4.1   Pore-scale seepage modeling based on LBM

Soil seepage processes are fundamentally governed by pore-scale fluid flow, interfacial tension, and structural connectivity, whereas macroscopic continuum models often rely on empirical constitutive relationships and struggle to reveal the intrinsic effects of heterogeneous pore structures on seepage behavior. As demonstrated by Martys and Chen (1996), the LBM can naturally recover equivalent permeability through statistical averaging of internal pore flows without explicitly introducing Darcy's law, thereby establishing a theoretical foundation for understanding soil seepage mechanisms from microscopic structures.

It has been common practice among many researchers to use one phase BGK or MRT LBM to systematically examine the flow properties of idealized and actual soil pore structures in case of saturated seepage condition. This was also discovered by Kang et al. (2002) and Patriarche et al. (2004); the pore size distribution width, connectivity of pores, and anisotropy are the main contributors to the effective permeability instead of porosity alone. The results of simulations of three-dimensional $X$-ray CT reconstructed soil pore structures indicate that even when it is identical porosity, the permeability of various soil structures can be several orders of magnitude apart, which contradicts the simple dominance of porosity in traditional macroscopic models as pointed out by Blunt et al. (2013), Lorenz et al. (2009) and Wang et al. (2023a). Such an outcome is of particular importance to comprehend the influence of tillage disturbance, compaction effects, and aggregate structure on infiltration rate in the field of hydro-mechanical behavior of soil.

Multiphase LBM models now capture water-air interface dynamics and capillary effects in unsaturated soils. The Shan-Chen pseudopotential model dominates this field. Liu et al. (2012) and Han et al. (2022) showed that pore-scale wettability patterns and local geometry jointly govern wetting front advance, producing fingering and preferential flow at larger scales. This finding departs from the Richards equation assumption of smooth wetting fronts. LBM simulations reveal pronounced front instability and strong dependence on initial moisture, particularly under rainfall or irrigation conditions.

The latest pore scale investigations also show that the concept of seepage behavior needs to be viewed as a structural and interfacial reaction, instead of a pore size dependent one. The information in Table 1 indicates that alterations in connectivity, phase distribution, and interface structure could significantly affect the local velocity field, permeability change and macroscopic asymmetry. Such mechanisms are especially apparent in the presence of unsaturated and freeze-thaw circumstances when interface movement and pore occlusion are constantly varying the actual seepage routes.

Table  1  Comparison of typical pore-scale soil seepage studies based on LBM

Study	Pore structure	Structure image	Flow type	Main conclusion
Shen et al. (2025)	QSGS + CT; 2 µm		Single-phase (with solute)	Water-salt migration is strongly coupled with pore heterogeneity; LBM is suitable for handling local concentration gradients
Xia et al. (2024)	Improved QSGS (icewater curved interface); 1 µm		Saturated single-phase (with phase change)	Freezing preferentially blocks main seepage channels, and permeability decay is dominated by connectivity degradation
Jiang et al. (2024)	Undisturbed CT; 150 µm		Unsaturated water-air two-phase	Pore connectivity and wall wettability jointly control wetting front propagation and fingering flow
Wang et al. (2023c)	SR-μCT; 2 µm		Unsteady single-phase	Pore-scale flow heterogeneity amplifies macroscopic anisotropy
Zhou et al. (2019)	SR-μCT; 3.7 µm		Saturated single-phase (with solute)	Biochar introduces bimodal pores, homogenizing pore-scale flow velocity while amplifying macroscopic dispersion
Zhang et al. (2016)	CT + overlapping spherical bubbles; 6 µm		Single-phase	Under the zero-shear-stress assumption at the water-air interface, discretization error of curved interfaces is the primary cause of microscale permeability underestimation

4.2   Water-heat-phase change and freeze-thaw structural evolution

The freeze-thaw process in porous media represents a complex interplay among phase change, heat transfer, and pore structure reconfiguration. Unlike isothermal seepage, freezing induces volume expansion of pore water (approximately 9%), generating internal stresses that alter pore connectivity while latent heat release modulates temperature gradients, creating a strongly coupled thermo-hydraulic system (Li and Yin, 2025). Figure 2 illustrates the pore-scale phase distribution during freezing, where ice preferentially nucleates in larger pores due to lower freezing point depression, effectively blocking primary drainage pathways while retaining unfrozen water films in micropores (Chuvilin et al., 2022; Watanabe and Flury, 2008).

Figure  2  Pore distribution characteristics of soil during freezing process (Li et al., 2024)

Download: Full-Size Img

On the mesoscale level, the lattice Boltzmann method can capture such phenomena using interconnected thermodynamic models. Enthalpy-based LBM uses latent heat as a source term of the thermal evolution equation and allows a direct simulation of moving phase interface without the need of any front tracking (Jiaung et al., 2001). Li et al. (2024) further developed this idea by adding the dynamic contact angle correction scheme on the pseudopotential model, which lowered the prediction error of ice-water interface curvature to 3.7%. An analysis of the results showed that there is a correlation between porosity and phase transformation kinetics: an increase in porosity, starting at 0.35 and ending at 0.48, led to a latent heat release rate enhancement of 22% indicating how pores control the movement speed of freezing fronts through thermally induced effects.

This additional aspect of solute exclusion during freezing makes the thermo-hydraulic response even more complicated. With the rejection of salt rejection during ice formation, localized concentration changes create osmotic potentials that cause extra moisture to move towards the frozen area-a process known as frost suction. The spatiotemporal variation of moisture and particle concentrations over time can be seen in Figure 3 showing how such a connection works (Wang et al., 2023b). They formulated a combined model involving dynamic permeability tensors which achieved good agreement with experimental data (R² = 0.91) between calculated moisture migration velocity and experimental data. This indicates that an increase in Cl⁻ concentration of 1 g/L will lower the phase change temperature by 0.28 K, thus forming a cycle whereby the buildup of salt leads to increased freezing point depression, enhanced movement of unfrozen water, and growth of ice lenses.

Figure  3  Spatiotemporal evolution of moisture and particle concentrations during soil pore-scale phase change process (Wang et al., 2023b)

Download: Full-Size Img

The shift from single freezing events to cyclic processes introduces structural memory effects that fundamentally alter soil behavior. Single freezing creates temporary ice blockage, whereas repeated freeze-thaw cycles induce irreversible pore structure reorganization through particle rearrangement, aggregate fragmentation, and crack formation (Leuther and Schlüter, 2021). Figure 4 illustrates these structural consequences. Xu et al. (2024) combined unidirectional freeze-thaw experiments with numerical simulations, revealing that cyclic freezing with self-weight loading produces silt enrichment zones at characteristic depths. These zones form because freeze-thaw cycling generates internal cracks that facilitate particle sorting, with finer particles migrating downward while coarser fractions remain stable—a process invisible to macroscopic continuum models but readily captured by pore-scale LBM simulation.

Figure  4  Effects of freeze-thaw cycling on soil particles under combined static loading and temperature conditions (Xu et al., 2024)

Download: Full-Size Img

Song et al. (2016) reconstructed soil structures using random generation-growth methods to simulate coupled heat-moisture transfer in unsaturated freezing soils, demonstrating that cyclic freezing produces hysteretic water retention curves distinct from primary drying or wetting. This hysteresis originates from ice-induced pore blockage that persists as metastable configurations through thawing, only partially recovering original connectivity. Klöffel et al. (2024) confirmed that post-thaw permeability often exceeds pre-freezing values due to crack formation, yet capillary characteristics degrade due to reduced pore connectivity—creating a decoupling between hydraulic conductivity and water retention that challenges traditional pedo-transfer relationships.

Current LBM frameworks face critical limitations in bridging these processes across scales. Static contact angle assumptions fail to capture dynamic interfacial tension variations during freezing, while thermo-mechanical coupling remains incomplete. Ji et al. (2021) demonstrated nonlinear stress-phase change rate relationships that remain unincorporated. Computational constraints persist: despite GPU acceleration achieving 0.18 ms per timestep (Ma et al., 2023), meter-scale simulations require 2.4×10⁵ grid cells, necessitating adaptive mesh refinement or hybrid LBM-DEM approaches (Younes et al., 2023). Future advances require integration of phase-field interface capturing for dynamic ice crystal growth, THMC coupling with explicit salt chemistry, and immersed boundary methods for evolving geometries (Qin et al., 2023). The convergence of LBM with μCT imaging (Liu et al., 2021; Zhang et al., 2025) offers particular promise for validating structure evolution predictions against real pore geometries, bridging the gap between mesoscopic mechanism and macroscopic engineering application.

4.3   Multiscale modeling and LBM upscaling issues

Multiscale modeling has become a central issue in LBM research on porous-media seepage and freeze-thaw processes. The main challenge is no longer resolving pore-scale phenomena themselves. It is determining how pore-scale results transfer reliably to larger scales. In this context, the key question is not how many scales can be listed. It is which structural and hydraulic indicators can be stably upscaled from pore-scale simulations to REV, Darcy-scale, or field-scale models.

CT-based digital cores provide an important basis for this effort. They preserve pore geometry, connectivity, and heterogeneity more realistically than idealized structures. Combined with LBM, they allow direct evaluation of permeability, tortuosity, and effective diffusivity. They also enable assessment of phase distribution and local flow heterogeneity under controlled pore-scale conditions (Fu et al., 2023; Sun et al., 2023; Wang et al., 2023a). However, current studies also show important limitations. The usefulness of such simulations depends strongly on image resolution. It is also sensitive to segmentation quality and REV representativeness. Furthermore, the robustness of computed parameters must be verified when the observation scale changes.

Andrade and Mital (2019) emphasized that soil processes span multiple orders of magnitude in scale. This characteristic makes single-scale models insufficient for capturing all relevant mechanisms. Existing scale-bridging strategies fall into three main categories. The first involves direct upscaling from pore-scale LBM simulations. These provide effective parameters for continuum models, such as permeability and hydraulic conductivity (Rentschler et al., 2022; Torquato, 2020). The second employs hybrid coupling with other numerical methods. DEM-LBM and LBM-FVM frameworks are particularly notable. They link pore-scale flow with particle rearrangement, deformation, or larger-scale heat and mass transport (Ma et al., 2024; Younes et al., 2023). The third uses data-driven or image-based multiscale reconstruction. In this approach, CT and digital core analysis combine with numerical simulation. The goal is to identify structural descriptors that bridge microscopic evolution and macroscopic behavior (Liu et al., 2023; Yin et al., 2023). Table 2 summarizes these modeling approaches, highlighting their respective advantages, limitations, and scales of application.

Although these advancements still exist, there are still great obstacles to overcome. Lourenço et al. (2024) studied use cases of LBMs when applied to the mass transfer process of multicomponent mixtures. This problem was recognized as one that cannot be predicted satisfactorily. The greatest issue is the interpretation of parameters throughout different scales. A micropore-based measure is not necessarily an efficient macroscopic descriptor. Many of the local variables used in seepage and frost-thaw processes are highly sensitive. Their reaction is strongly influenced by surface properties and transient connection failure. These measures are sensitive to changes due to the presence of salinity as well as path-dependent structural rearrangements (Han et al., 2024; Li et al., 2024). As a result, parameter upscaling usually turns out to be unreliable. Predictive models might rely entirely upon chosen REV. They may also change with reconstruction method or coupling strategy.

The future work needs to be redirected. Rather than just expanding the size of simulation domains, it should direct its efforts towards finding indicators at the pore scale that can maintain a physical significance during upscaled stages. Such an effort will have particular relevance in respect of frozen and coastal porous formations where both the boundary conditions and the phase behavior tend to be more complicated. The objective can only be reached through better linkage between LBM and quantitative imaging technologies. It further depends on strict evaluation based on experiments at different scales (Andrade and Mital, 2019; Lourenço et al., 2024).

5   Conclusion

The LBM provides a mesoscale framework for studying seepage and freeze-thaw processes. Its advantage is not limited to spatial resolution. It can represent pore geometry, multiphase interfaces, latent heat effects, and connectivity changes within a unified model. Current studies show that the hydraulic response of frozen soil is not controlled by pore size alone. It is more strongly influenced by pore connectivity, interface dynamics, and unfrozen water migration. This helps explain why pore-scale LBM can capture mechanisms that are not resolved by continuum models. However, several limitations still restrict its application. These include limited computational scale, unclear physical interpretation of model parameters, difficulties in transferring pore-scale results to larger scales, and a lack of experimental validation. Future work should focus on improving the description of dynamic interfaces and phase-change processes. More reliable methods for scale bridging and parameter upscaling are also needed. In addition, coupled THMC models require further development, supported by imaging data and experimental observations. These challenges are more pronounced in coastal and marine-influenced porous media. In such environments, salinity changes, interactions among sea ice, permafrost, and groundwater, and tidal effects all influence freezing temperature, interface behavior, and seepage pathways. As a result, conclusions derived from inland frozen soils cannot be directly applied to these systems.