Multi-Phase Fluid Transport in Expansive Soils: A Fractal-Constrained Physics-Informed Neural Network Framework for Permeability and Capillary Flow Characterization

1 Introduction

Expansive soils are widely distributed worldwide and frequently host rainfall-induced shallow landslides, posing persistent threats to highways, embankments, and natural slopes [1,2]. The mechanism underlying these failures is fundamentally a problem of two-phase (water–air) fluid transport through unsaturated, micro-structured pore networks under transient hydraulic loading [3]. The water-retention behavior of compacted expansive clays, in particular, is highly sensitive to the dual-porosity micro-structure formed during compaction [4,5], and its evolution under wetting–drying cycles further complicates seepage prediction [2,6]. Stabilization measures such as additive treatment have been studied to mitigate desiccation and cracking [7,8], yet the unsaturated hydraulic response remains the most uncertain input in slope stability analysis.

As a core factor affecting the accuracy of slope stability evaluation, soil parameter variability has been widely concerned in probabilistic slope analysis, Jiang et al. [9]. Constitutive modelling of unsaturated, state-dependent soil behavior has matured considerably over the past two decades [10,11], leading to widely adopted frameworks such as the van Genuchten (VG) model. In practice, however, the VG parameters α and n are usually extracted by curve-fitting the capillary pressure–saturation (Pc–Sw) data obtained from mercury intrusion porosimetry (MIP). This procedure carries a methodological blind spot: MIP is itself a fluid displacement experiment governed by the Laplace–Washburn balance. Mercury, as a non-wetting fluid, invades the pore network under controlled pressure in a process that is mechanically analogous to water displacing air during rainfall infiltration. Treating MIP solely as a structural probe—while ignoring its fluid-displacement nature—represents a methodological missed opportunity. The present study attempts to reverse that perspective: MIP data are reinterpreted as records of a multi-phase fluid displacement experiment, from which the parameters that actually govern fluid transport are inverted, rather than merely fitting a curve.

A second under-exploited aspect of MIP data lies in the fractal description of the pore architecture. While existing fractal analyses [12,13,14,15] mainly treat the fractal dimension Ds as a geometric descriptor of surface roughness, from a fluid-mechanics standpoint, Ds is more meaningfully understood as an indicator of pore surface complexity and structural heterogeneity, which generally implies more complex potential flow paths. It is important to clarify, however, that Ds is not a direct measure of hydraulic tortuosity. The classical relationship τ = ( L 0 2 r ) ( D f − 1 ) relates tortuosity to the pore volume fractal dimension Df, whereas our MIP data primarily yield the surface fractal dimension Ds. While Ds and Df are theoretically related in some frameworks, the connection is not mathematically exact for all pore geometries. We therefore interpret Ds as a robust relative indicator of structural complexity that correlates with, but does not strictly equate to, hydraulic tortuosity. This correlation provides the theoretical justification for using the trend in Ds as a physics-informed prior constraint in the FC-PINN inversion. Four established fractal models are independently applied to ensure robustness: the Zhang–Li model [12]; the Xu model [13]; the Friesen–Mikula model [14]; and the Neimark model [15].

Physics-informed neural networks (PINNs) are particularly well suited to performing this bridging task. Unlike purely data-driven methods, PINNs embed governing equations into the loss function, requiring the network to satisfy both data fidelity and physical constraints simultaneously [16]. Subsequent reviews [17,18] have established PINNs as a powerful paradigm for inverting partial differential equations of transport phenomena. In the context of soil-water flow, Bandai and Ghezzehei [19] developed monotonicity-constrained PINNs to estimate constitutive relationships and flux density from infiltration data; He et al. [20] coupled PINNs with multiphysics data assimilation for subsurface transport; and Wang et al. [21] proposed theory-guided neural networks for subsurface flow modelling. More recently, Depina et al. [22] applied PINNs to inverse problems in unsaturated groundwater flow, while Shadab et al. [23] employed PINNs to investigate steady unconfined groundwater flow and invert hydraulic conductivity from hydrologic observations. Ng et al. [24] developed a novel surrogate model for hydro-mechanical coupling in unsaturated soil with incomplete physical constraints, demonstrating the growing maturity of PINN-based approaches in geotechnical engineering. Most recently, Zhao et al. [25] demonstrated PINN-based numerical simulation of multi-phase flow mass transfer, signalling the maturity of PINNs for multi-phase problems. To the authors’ knowledge, however, no prior study has applied PINNs to multi-phase flow inversion in expansive soils, nor has any prior work incorporated MIP-derived fractal dimensions as a physics-informed prior within the PINN loss structure.

The contributions of this work are fourfold: (i) MIP data are reinterpreted from a fluid-mechanics perspective, establishing their mechanical analogy with unsaturated water–air flow; (ii) a multi-objective FC-PINN loss function is constructed in which the steady-state Richards equation residual is the core physics term, regularized by fractal-dimension priors derived from four benchmark MIP models [12,13,14,15]; (iii) the framework is systematically validated against three void-ratio specimens covering the practical compaction range; (iv) the topological mechanism behind the steep dependence of permeability on void ratio is elucidated, with implications for slope seepage analysis and rainfall-triggered landslide early warning.

2 MIP as a Fluid Displacement Experiment and Fractal Analysis

2.1 Test Material

The expansive soil used in this study was sampled from the Quaternary Jingmen Formation in Hubei Province, China, and classified as high-plasticity clay according to the Unified Soil Classification System. Index properties include a liquid limit of 67.4%, a plastic limit of 29.8% (plasticity index 37.6%), a specific gravity Gs = 2.74, and a free-swell ratio of 68%—consistent with the indices reported for expansive soils in adjacent regions [1,6]. Cylindrical specimens (Φ35 × 70 mm) were prepared by static layered compaction to target void ratios e = 1.4, 1.5 and 1.6 (designated S1, S2, S3) and equilibrated for 48 h in a humidity-controlled chamber prior to MIP testing. A Micromeritics AutoPore IV 9520 porosimeter was used, with a pressure range of 0.003–227 MPa corresponding to pore-throat diameters from 360 μm down to 5 nm.

2.2 Mechanical Isomorphism of Fluid Displacement

In MIP, mercury invades the pore network as a non-wetting phase, with the driving pressure P related to the throat radius r by the Washburn equation: P=−2σcosθcr(1) where σ = 0.485 N/m is the surface tension of mercury and θc = 130° is the mercury–soil contact angle (note that cos130° is negative, hence the negative sign ensures positive pressure). This functional form is identical to the capillary pressure–throat radius relationship for water–air systems in unsaturated soils [4,5]. From a fluid-mechanics viewpoint, MIP is therefore an accelerated analog of water–air displacement: although the fluids and time-scales differ, the underlying force balance—the competition between capillary forces and the externally applied pressure differential—is conceptually analogous. This analogy is the physical foundation underpinning every subsequent argument in this paper (Fig. 1).

It is crucial to acknowledge the inherent limitations of this analogy. First, differences in wettability (mercury as a non-wetting phase versus water as a wetting phase) lead to different micro-scale displacement mechanisms. Second, contact angle hysteresis and “ink-bottle” effects mean that MIP-derived pore size distributions are often biased towards pore throat sizes rather than true pore body sizes. Third, the high-pressure intrusion process in MIP operates on vastly different time scales and mechanical paths compared to low-velocity, transient natural rainfall infiltration. Therefore, our framework treats MIP as an accelerated mechanical analog whose value lies in providing prior constraints based on thermodynamic and fluid mechanical similarities, rather than establishing a point-by-point physical equivalence. This careful positioning ensures that the subsequent FC-PINN inversion is understood as an analogy-driven, data-informed strategy.

Figure 1: Fluid-mechanics-based reinterpretation of MIP: mercury intrusion as an accelerated mechanical analog of water–air displacement in unsaturated soils. The analogy is grounded in the common Laplace–Washburn force balance, though differences in wettability, contact angle hysteresis, ink-bottle effects, and loading rates must be acknowledged as inherent limitations.

2.3 Fractal Dimension as a Structural Complexity Indicator

Conventional analyses treat the fractal dimension Ds as a geometric parameter [12,13,14,15]. From a fluid-mechanics standpoint, however, Ds is more meaningfully understood as an indicator of pore surface complexity and structural heterogeneity. It is important to clarify that surface fractal dimension and hydraulic tortuosity are distinct concepts. The classical relationship τ = ( L 0 2 r ) ( D f − 1 ) relates tortuosity to the pore volume fractal dimension Df, whereas our MIP data primarily yield the surface fractal dimension Ds. While Ds and Df are theoretically related (e.g., Ds = Df + 1 in some frameworks), the connection is not mathematically exact for all pore geometries. We therefore interpret Ds as a robust relative indicator of structural complexity that correlates with, but does not strictly equate to, hydraulic tortuosity. This correlation provides the theoretical justification for using the trend in Ds as a physics-informed prior constraint in the FC-PINN inversion, as detailed in Section 3.2. Four established fractal models are independently applied to ensure robustness: (i) the Zhang–Li model (regression of cumulative intrusion work [12]); (ii) the Xu model (cumulative pore volume versus radius [13]); (iii) the Friesen–Mikula model [14]; and (iv) the Neimark model based on surface-energy integration [15].

2.4 Conversion of Mercury-Intrusion Data to P_c–S_W Dataset

Capillary pressure Pc for each measurement point is directly calculated by the Washburn equation (Eq. (1)) from applied intrusion pressure and pore-throat radius.

The cumulative intruded mercury volume VHg(P) at pressure P is converted to water saturation Sw using the total sample pore volume Vp derived from specimen bulk density and specific-gravity: V p = m s ( 1 ρ b − 1 G s ρ w ) , where ms is solid mass of specimen; ρ b is bulk density; Gs is specific gravity of soil solids; ρ w is water density.

Mercury porosimetry cannot access intra-aggregate nanoscale pores below the instrument minimum pore-size cut-off (5 nm). This inaccessible pore fraction is treated as permanently water-filled (residual saturation Sw,r. The accessible pore volume Vp,acc is the maximum measured cumulative mercury-intrusion volume: V p , a c c = m a x [ V H g ( P ) ] .

The degree of non-wetting-phase mercury saturation SHg is: S H g ( P ) = V H g ( P ) V p , a c c .

By the pore-space occupancy analogy for non-wetting/wetting phase, the equivalent water saturation for the same pore-throat configuration is: S w ( P ) = S w , r + ( 1 − S w , r ) · ( 1 − S H g ( P ) ) .

Trapped mercury after pressure release is not used in this calculation: only in-process cumulative intrusion volume under applied pressure is adopted, consistent with the MIP pressure-volume curve recorded during loading.

Ink-bottle effects systematically bias measured pore-throat radii towards smaller values; we do not perform explicit ink-bottle correction. Instead, this structural bias is treated as inherent noise of the experimental dataset, and the PINN’s physical regularisation (Richards-equation residual + fractal prior) acts to regularise and smooth these noisy Pc–Sw inputs during inversion. The fitted VG parameters therefore represent effective, ink-bottle-affected equivalent parameters for the tested compacted expansive soil specimens, not ideal intrinsic pore-geometry constants.

3 FC-PINN Framework: Embedding Fluid Equations into Learning

3.1 Governing Equations

To avoid ambiguity, we first define all variables and their relationships. The primary variables are:

P [Pa]: Mercury intrusion pressure in MIP.

Pc [Pa]: Capillary pressure, defined as Pc = Pair − Pwater. In this study, Pc is derived from MIP data via the Washburn equation (Eq. (1)).

ψ [m]: Pressure head, related to capillary pressure by ψ = − P c ρ w · g , where ρw is the density of water and g is gravitational acceleration.

Sw [-]: Water saturation, the fraction of pore volume occupied by water.

Se [-]: Effective saturation, defined as S e = S W − S w , r S w , s − S w , r , where Sw,r and Sw,s are the residual and saturated water saturations, respectively.

θ [-]: Volumetric water content, θ = n·Sw, where n is the porosity. Residual and saturated water contents are θr = n·Sw,r and θs = n·Sw,s.

kr [-]: Relative permeability, a function of Se.

K(Se) [m/s]: Unsaturated hydraulic conductivity, K(Se) = Ks·kr(Se), where Ks is the saturated hydraulic conductivity.

Unsaturated water–air flow in expansive soil slopes is governed by Richards’ equation [19,20,21]:

In Eq. (2), θ is the volumetric water content [-]; K(Se) is the unsaturated hydraulic conductivity [LT−1] (e.g., m/s), which is a function of the effective saturation Se; Ψ is the pressure head [L]; ê z is the unit vector in the vertical direction (with the z-coordinate taken as positive upward); and qs is the source/sink term [T−1] (here qs = 0 for the closed system). The VG constitutive model [26] is adopted to describe the soil hydraulic properties [10,11].

In this work, the van-Genuchten parameter α is expressed in MPa−1 to match the pressure units of MIP-derived capillary-pressure data, rather than the conventional pressure-head unit m−1.

It is important to clarify the role of the Richards equation in the present framework. Unlike conventional PINN applications that solve transient Partial Differential Equations with spatial and temporal coordinates as inputs, our FC-PINN is designed as a physically regularized neural-network curve-fitting approach for inverting steady-state or quasi-static constitutive relationships from MIP data. The physical residual is therefore evaluated using the steady-state form of Eq. (2): ∇ · K S e ∇ Ψ + ê z = 0 . This residual is computed in a transformed coordinate system where log10(Pc) serves as a proxy for spatial location, ensuring the inferred Se(ψ) and K(Se) relationships are physically consistent.

3.2 Architecture and Loss Function

Following the original PINN formulation of Raissi et al. [16] and the subsequent multiphysics extensions [17,18,20], a fully connected network with 5 hidden layers of 64 neurons each is constructed, using tanh activation functions. Inputs are the normalized log10(Pc) and the void ratio e; outputs are Sw and the relative permeability kr. The VG parameters α, n and Ks are treated as trainable physical parameters, reparameterized through sigmoid functions to enforce physically reasonable bounds. To mitigate overfitting given the modest dataset size, we employ early stopping (training ceases when the validation loss does not decrease for 200 consecutive epochs) and L2 weight regularization. Sensitivity analysis on network architecture (testing 3 × 32, 4 × 48, and 5 × 64 configurations) indicated that 5 × 64 achieved the best balance between accuracy and generalization, with a training–validation loss gap consistently below 8%. The total loss is:

All loss-component terms are computed on network-normalised variables; raw physical parameters are recovered via inverse scaling for tabulation and interpretation.

Where the individual loss terms are defined as follows:

Ldata: The mean squared error between the FC-PINN predictions and the MIP-derived Sw data. This is the data-fidelity term.

Lphys: The physics residual. It enforces the steady-state Richards equation. To compute this residual, we use a 1D vertical coordinate transformation. Recognizing that the capillary pressure head ψ is a natural spatial variable in unsaturated flow, we introduce a transformed surrogate coordinate x = log 10 ( − ψ ) = log 10 ( P c ρ w · g ) . which is derived from pressure rather than a directly spatial coordinate. This transformation is justified under the one-dimensional hydrostatic assumption, where the gradient in log 10 ( P c ) directly corresponds to a spatial gradient in elevation. The residual is then discretized using automatic differentiation to compute the derivatives of Se and K with respect to x. The explicit expression for the loss is: Lphys=1Nx∑i=1Nx(ddxK(Se(xi))(dψdx+1))2(6) here, Nx is the number of collocation points in the transformed spatial domain. This approach avoids the need for explicit spatial coordinates as network inputs.

Lfractal: The fractal prior loss. Based on the empirical correlation observed between the fractal dimension Ds and the VG parameter n (i.e., higher Ds corresponds to a broader pore-size distribution and thus a lower n), we construct a weakly-constrained loss term. Instead of a hard constraint, we define Lfractal as a regularization term that penalizes deviations of the inferred n from the trend established by the Ds values. For each specimen, the target ntarget is estimated via linear regression between the mean Ds from the four models (Table 1) and the n values from a preliminary conventional VG curve fit (denoted as ninitial). The regression equation is n t a r g e t = a · D s ¯ + b . With a = −3.76 and b = 11.84 (coefficient of determination R2 = 0.97). The individual mean D s ¯ and corresponding preliminary ninitialninitial values used to determine this regression are provided in Table 2. Lfractal is then defined as:

For numerical stability inside the PINN, both n and n t a r g e t are linearly mapped into a dimensionless normalised interval [0, 1] before computing L f r a c t a l = ( n − n t a r g e t ) 2 . All reported n values in tables are converted back to raw physical units after inversion. The loss weights λ1–λ4 are tuned for this normalised objective space.

D s ¯ is the arithmetic mean of the four model-specific Ds values listed in Table 1 of the main text (Zhang–Li, Xu, Friesen–Mikula, and Neimark models).

ninitial is obtained from a preliminary conventional VG curve fit to the raw MIP data, serving as the dependent variable for establishing the empirical fractal–n trend.

ntarget is computed using the regression equation ntarget = −3.76· D s ¯ + 11.84, which serves as the weakly constrained physical prior in the fractal loss term Lfractal = (n − ntarget)2.

This approach is consistent with the interpretation of Ds as a relative indicator of structural complexity, and guides the PINN inversion towards physically consistent parameter space without imposing an overly rigid constraint.

Lbc: The boundary condition loss. It enforces physical asymptotes: as Pc → 0 (i.e., saturation), Sw → Sw,s; and as Pc → ∞ (i.e., residual), Sw → Sw,r. The boundary points are selected as the minimum and maximum Pc values from the MIP dataset. The loss is: Lbc=(Sw(Pc,min)−Sw,s)2+(Sw(Pc,max)−Sw,r)2(8) where Sw,s and Sw,r are determined from the MIP data asymptotes as described in Section 4.1.

The weights λ1 = 1.0, λ2 = 0.5, λ3 = 0.3, λ4 = 0.2 were obtained from a grid search on a validation set, consistent with the weighting strategies suggested in recent PINN reviews [17,18]. Sensitivity analysis showed that the inferred VG parameters varied by less than 5% when weights were perturbed by ±20%, confirming robustness to weight selection. Training proceeded in two stages: 2000 Adam iterations for rapid convergence, followed by L-BFGS-B for fine refinement. The initial learning rate was 10−3, decayed by a factor of 0.5 every 500 epochs. For each specimen, the MIP-derived Pc–Sw dataset (approximately 150 points) was split via stratified random sampling into 80% training (≈120 points) and 20% validation (≈30 points) The architecture and loss structure are illustrated in Fig. 2.

Figure 2: Architecture of the Fractal-Constrained Physics-Informed Neural Network (FC-PINN) framework with composite loss function combining MIP data fidelity, steady-state Richards equation residual, and fractal-prior regularization. The network inputs are log₁₀(P_c) and void ratio e; outputs are S_w and K_r. The VG parameters α, n, and K_s are trainable physical parameters embedded within the loss.

4 Results and Discussion

4.1 Pore-Size Distribution and Capillary Pressure–Saturation Curves

The MIP-derived pore-size distributions of S1, S2 and S3 are shown in Fig. 3a. All specimens exhibit characteristic bimodal distributions: a macropore peak at 15–20 μm (inter-aggregate pores) and a micropore peak at 0.05–0.1 μm (intra-aggregate pores). This pattern is consistent with the dual-porosity structure widely reported for compacted clayey soils. As the void ratio rises from 1.4 to 1.6, the macropore peak grows significantly (peak incremental intrusion increases from 0.0116 to 0.0148 mL/g), whereas the micropore peak remains essentially unchanged—indicating that compaction state predominantly affects the inter-aggregate porosity. The total intrusion volumes are 0.224, 0.261 and 0.298 mL/g, corresponding to MIP-estimated porosities of 0.381, 0.417 and 0.448, respectively. From these data, the residual and saturated water saturations (Sw,r and Sw,s) are determined as follows: Sw,s is taken as the maximum Sw value at the lowest Pc (approaching 1.0), and Sw,r is taken as the minimum Sw value at the highest Pc (approaching 0.0). The effective saturation Se is then calculated using the relationship defined in Section 3.1.

Fig. 3b presents the Pc–Sw curves derived from Neimark’s surface-energy method [15] together with the FC-PINN fit. The air-entry pressures are approximately 0.85, 0.71 and 0.58 MPa for S1, S2 and S3, decreasing systematically with void ratio. We highlight the intermediate saturation range Sw = 0.3–0.7, which corresponds to the most hydraulically active regime during rainfall infiltration into slopes—the predictive accuracy in this band directly determines the reliability of any landslide early warning system.

Figure 3: (a) MIP pore-size distributions for three void ratios (S1: e = 1.4, S2: e = 1.5, S3: e = 1.6); (b) capillary pressure–saturation curves with FC-PINN fits. The highlighted band (S_w = 0.3–0.7) marks the hydraulically active saturation regime critical for rainfall infiltration analysis.

4.2 Fractal Dimension and Flow-Path Complexity

The Ds values from the four models are summarized in Table 1, with the Zhang–Li log–log regression shown in Fig. 4a. All regressions exhibit excellent linearity (the mean R2 values range from 0.9944 to 0.9981), confirming the self-similar fractal character of the pore surfaces over the measurement range [14]. Ds increases monotonically with void ratio. This ~3% change in Ds appears modest, yet the corresponding permeability changes by nearly fourfold—revealing that the dramatic permeability response is not driven by mean pore size, but by a fundamental reorganization of the pore connectivity topology [27].

As discussed in Section 2.3, Ds is interpreted here as an indicator of structural complexity. The systematic increase in Ds with void ratio reflects greater surface irregularity and pore-structure heterogeneity. However, the modest ~3% increase in Ds contrasts sharply with the ~4× increase in Ks, suggesting that permeability is governed primarily by topological connectivity—specifically, the increased volume and connectivity of inter-aggregate macropores—rather than by pore-scale roughness or tortuosity alone.

Fig. 4b presents the Ds values from all four models alongside the FC-PINN-inverted saturated permeability Ks as functions of void ratio. While the absolute Ds values differ systematically across models (Neimark consistently yields the highest values), the trend with void ratio is fully consistent. The relative variation in fractal dimension is thus a robust indicator of structural complexity and pore-network heterogeneity, irrespective of which model is adopted.

Figure 4: (a) Zhang–Li log–log fractal regression for three specimens; (b) fractal dimensions from four MIP-based models and PINN-inverted saturated permeability versus void ratio.

4.3 FC-PINN Training Performance and Parameter Inversion

The FC-PINN convergence trajectory is shown in Fig. 5a. The Richards-equation residual Lphys decreases rapidly from ~0.75 to ~0.05 within the first 700 epochs, reflecting efficient learning of the fluid physics. The data-fidelity loss Ldata plateaus at ~0.02 after about 1500 epochs. The L-BFGS-B refinement stage further reduces the total loss by approximately 23%. Final total losses for S1, S2, S3 are 0.0041, 0.0037 and 0.0033, respectively, and the training–validation loss gap remains below 8%, indicating no overfitting. The dataset split (80/20 training/validation) for each specimen was performed via stratified random sampling, ensuring that both sets cover the full Pc–Sw range. The small and consistent training–validation gap confirms that the chosen architecture (5 × 64) with early stopping and L2 regularization effectively prevents overfitting despite the modest dataset size. Compared with the monotonicity-constrained PINN approach of Bandai and Ghezzehei [19], the fractal-prior constraint adopted here produces tighter physical bounds on the inferred VG parameters without sacrificing fit quality.

The FC-PINN-inverted VG parameters are listed in Table 3 and admit clear fluid-physical interpretations: (i) α increases from 1.42 to 2.18 m−1, signifying a lower air-entry suction and easier destabilization of the water–air interface in larger pores—consistent with the expansion of inter-aggregate porosity at higher void ratios [4,5]; (ii) n decreases from 1.67 to 1.51, indicating a broader pore-size distribution, so that fluids encounter both fast flow channels and stagnant zones simultaneously, potentially triggering non-Darcian effects and local non-equilibrium flow; (iii) Ks rises by nearly a factor of four while Ds increases by only ~3%, a striking nonlinearity that reveals the topological—not metric—origin of the permeability sensitivity to void ratio. Specifically, this suggests that the primary mechanism for the permeability increase is the creation of additional inter-aggregate macropore pathways (topological connectivity) rather than a reduction in flow-path tortuosity per se.

FC-PINN-retrieved VG hydraulic parameters are compiled in Table 3. With α given in MPa−1, the characteristic air-entry suction 1/α ranges from 0.46 to 0.70 MPa, which is physically consistent with the measured air-entry pressures (0.58–0.85 MPa) obtained from MIP Pc-Sw curves. Minor deviations arise from ink-bottle effects and the soft fractal regularisation during PINN inversion. Clear physical interpretation can be drawn: (i) α increases from 1.42 MPa−1 to 2.18 MPa−1, corresponding to reduced air-entry suction and easier water-air interface destabilisation within larger inter-aggregate pores at higher void ratios [4,5]; (ii) n decreases from 1.67 down to 1.51, indicating broader pore-size distribution; co-existence of large fast-flow channels and stagnant small pores potentially induces non-Darcian behaviour and local flow non-equilibrium; (iii) Ks rises almost four-fold while Ds only increases by ~3%. This strong nonlinearity proves that void-ratio-driven permeability amplification originates from topological connectivity changes (creation of new inter-aggregate macropore flow pathways), rather than simple reduction of flow-path tortuosity.

Figure 5: (a) FC-PINN training-loss convergence for sample S1; (b) RMSE comparison of four models across three specimens. The fractal-constrained PINN shows superior accuracy and lower variability. Zhang–Li [12], Xu [13], Friesen–Mikula [14].

4.4 Comparison with Conventional Models

Table 4 and Fig. 5b quantitatively compare the FC-PINN with three conventional curve-fitting models [12,13,14]. The FC-PINN achieves a mean root mean square error (RMSE) of 0.0198 and a mean R2 of 0.9952. Relative to the best conventional model (Xu), this represents a 31% reduction in RMSE and a 0.85% improvement in R2. Of greater engineering relevance: in the intermediate saturation range (Sw = 0.3–0.7), where the fluid is most active during rainfall infiltration, the PINN reduces the error by 42%. The coefficient of variation of RMSE across the three specimens is only 3.6% for the FC-PINN versus 8.1–14.2% for the conventional models, demonstrating substantially greater robustness—an observation consistent with the multi-phase PINN performance reported by Zhao et al. [25].

4.5 Implications for Slope Seepage Modelling

The Sw–Pc and kr–Se relationships from the FC-PINN can be directly imported into Richards-equation-based slope seepage solvers to perform simulations that are far more physically faithful than conventional fits would allow. Fig. 6 presents a conceptual slope model with heterogeneous void-ratio zones, motivated by the spatially variable soil-property analyses of Jiang et al. [9], together with the relative permeability curves predicted using the FC-PINN-inverted n values. Two engineering consequences emerge (Fig. 6a): (1) in high-e zones (e.g., the loose toe region, e = 1.6), n is small → kr rises steeply near saturation (Fig. 6b) → rainfall quickly drives the soil toward full saturation → matric suction drops abruptly → slope stability deteriorates rapidly, raising the risk of landslide triggering; (2) in low-e zones (e.g., the compacted crest, e = 1.4), n is larger → fluid migration is slower → transient perched water tables can form, locally elevating pore-water pressure as an alternative triggering mechanism. These observations align with the three-dimensional crack-propagation patterns reported by Zhang et al. [2] for expansive soil slopes under wetting cycles, and provide a physically interpretable basis for the soilbag-protected slope design promoted by Xu and Zhang [1].

Figure 6: (a) Conceptual slope profile with heterogeneous void-ratio zones (crest: e = 1.4, mid-slope: e = 1.5, toe: e = 1.6); (b) Mualem k_r–S_e curves derived from FC-PINN-inverted n values for the three compaction states. The curves illustrate how the n parameter governs the steepness of the relative permeability response near saturation, with direct implications for slope seepage behavior.

5 Conclusions

This study reinterprets MIP data on expansive soils from a fluid-mechanics perspective and develops a Richards-equation-constrained, fractal-prior-informed PINN framework for multi-phase flow parameter inversion. The principal conclusions are as follows:

• (1)MIP is, by its very mechanics, a multi-phase fluid displacement experiment that is isomorphic to water–air flow in unsaturated soils. Making this isomorphism explicit transforms MIP from a passive structural probe into an active record of fluid displacement [4,5].
• (2)The three expansive soil specimens (e = 1.4, 1.5, 1.6) exhibit bimodal pore-size distributions with macropore peaks at 15–20 μm and micropore peaks at 0.05–0.1 μm. Inter-aggregate porosity expands substantially with void ratio, and air-entry pressure decreases from 0.85 to 0.58 MPa.
• (3)The four fractal models [12,13,14,15] yield Ds values that are systematically correlated and increase monotonically (2.47 → 2.55). The robust trend in relative Ds supports its use as a quantitative indicator of structural complexity that correlates with flow-path tortuosity.
• (4)The fractal-constrained PINN simultaneously satisfies data fidelity, the Richards equation residual, and the fractal prior, returning physically meaningful α, n, Ks values. Compared with the best conventional model, RMSE is reduced by 31% and R2 reaches 0.9952; in the engineering-critical intermediate saturation range, the error drops by 42%.
• (5)The four-fold increase in permeability accompanied by only ~3% increase in Ds reveals that permeability sensitivity to void ratio is rooted in topological reorganization of pore connectivity rather than in mean pore size. This provides a new, physically interpretable indicator for slope-stability prediction and rainfall-triggered landslide early warning [2,9].