Physics-constrained symbolic regression for discovering closed-form equations of multimodal water retention curves from experimental data

Yejin Kim; Hyoung Suk Suh

physics.geo-phFeb 24, 2026

Physics-constrained symbolic regression for discovering closed-form equations of multimodal water retention curves from experimental data

Authors: Yejin Kim, Hyoung Suk Suh

Organizations: Department of Civil and Environmental Engineering, Case Western Reserve University, United States

Abstract

Modeling the unsaturated behavior of porous materials with multimodal pore size distributions presents significant challenges, as standard hydraulic models often fail to capture their complex, multi-scale characteristics. A common workaround involves superposing unimodal retention functions, each tailored to a specific pore size range; however, this approach requires separate parameter identification for each mode, which limits interpretability and generalizability, especially in data-sparse scenarios. In this work, we introduce a fundamentally different approach: a physics-constrained machine learning framework designed for meta-modeling, enabling the automatic discovery of closed-form mathematical expressions for multimodal water retention curves directly from experimental data. Mathematical expressions are represented as binary trees and evolved via genetic programming, while physical constraints are embedded into the loss function to guide the symbolic regressor toward solutions that are physically consistent and mathematically robust. Our results demonstrate that the proposed framework can discover closed-form equations that effectively represent the water retention characteristics of porous materials with varying pore structures. To support third-party validation, application, and extension, we make the full implementation publicly available in an open-source repository.

Figures & tables

Figure 1: Schematic illustration of the data mapping procedure. Black hollow symbols represent exemplary water retention data, with the black axes denoting the original

(s,S_{w})

space and the red axes indicating the mapped

(s^{*},S^{*}_{w})

space. The upper left corner of the red box corresponds to the reference point

(s_{\text{min}},S_{w,\text{max}})

and the lower right corner to

(s_{\text{res}},S_{w,\text{min}})

in the original space, which map to

(0,1)

and

(1,0)

, respectively, in the transformed space. This mapping ensures that the learned function (red curve) resides within

s^{*}\in[0,1]

and

S^{*}_{w}\in[0,1]

Figure 2: Exemplary symbolic equation represented as a binary tree. This tree represents an expression $(1+x)\ln(x)-\exp(y)+0.5$ , with a depth of 4 and a total size of 12 nodes.

Panel of Figure 3 — Figure 3: Illustration of (a) mutation and (b) crossover operations in an evolutionary symbolic regression algorithm.

Panel of Figure 4 — Figure 4: The symbolic model discovered from experimental unimodal water retention data using the proposed approach, compared against the ( Van Genuchten, 1980 ) model, vanilla SR, and PCSR without mode constraint.

Panel of Figure 5 — Figure 5: Data and physics loss components and the number of modes exhibited by PCSR with $\mathcal{L}_{\text{mode}}$ at varying complexity levels that optimally fit the data presented in Figure 4 .

Panel of Figure 6 — Figure 6: Learning curves of loss components and predicted number of modes for the PCSR without and with $\mathcal{L}_{\text{mode}}$ for the data presented in Figure 4 .

Panel of Figure 7 — Figure 7: The symbolic model discovered from experimental bimodal water retention data using the proposed approach, compared against the ( Durner, 1994 ) model, vanilla SR, and PCSR without mode constraint.

Panel of Figure 8 — Figure 8: The symbolic model discovered from synthetically generated multimodal water retention data, compared against the ( Durner, 1994 ) model, vanilla SR, and PCSR without mode constraint.

Appendix figures & tables6 assets

Supplementary material from the paper’s appendix.

Appendix

Model	Discovered expression
Vanilla SR	$\begin{aligned} y&=2\times x\times(1.5181\times\exp(x\times(-0.0948\times x+6.829\times 10^{-6}))\\ &+(0.0327\times x\times x-9.953\times 10^{-6}\times x-0.4680))+1.5125\times 10^{-6}\end{aligned}$
SR with monotonicity cstr.	$\begin{aligned} y&=((\exp(((\exp(((x\times(-0.176))\times(x+(\exp((x\times x)\times 0.469)\\ &\times(x\times 0.434))))+0.465)+((x+x)\times x))\\ &\times(-0.108))+0.502)+0.718)\times x)-0.001\end{aligned}$

Appendix

Table A.1 : Discovered equations that best fit the data generated from the benchmark expression in Eq. ( A.1 ), corresponding to those illustrated in Figure A.1 .

Panel of Figure 10 — Figure A.1 : Effect of the monotonicity constraint ( ${\mathrm{d}y}/{\mathrm{d}x}\geq 0$ ) in the noise-free case. (a) Benchmark equation and sampled data (black hollow circles) compared with vanilla SR (green solid curves) and SR with monotonicity constraint (red solid curves). (b) First derivatives of the corresponding equations, illustrating monotonicity violations in vanilla SR and preservation in the constrained model.

Panel of Figure 11 — Figure A.2 : Effect of the monotonicity constraint ( ${\mathrm{d}y}/{\mathrm{d}x}\geq 0$ ) in the presence of Gaussian process noise. (a) Noisy sampled data and learned equations from vanilla SR and SR with monotonicity constraint. (b) First derivatives of the corresponding equations, illustrating monotonicity violations in vanilla SR and preservation in the constrained model.

Panel of Figure 12 — Figure A.3 : Effect of the limiting constraint ( $y=0$ for $x<0.2$ ) in the noise-free case. (a) Benchmark equation and sampled data compared with vanilla SR and SR with limiting constraint. (b) First derivatives of the corresponding equations, illustrating strict enforcement of the limiting behavior in the constrained model.

Panel of Figure 13 — Figure A.4 : Effect of the limiting constraint ( $y=0$ for $x<0.2$ ) in the presence of Gaussian process noise. (a) Noisy sampled data and learned equations from vanilla SR and SR with limiting constraint. (b) First derivatives of the corresponding equations, demonstrating reduced sensitivity to noise under the imposed limiting constraint.

Panel of Figure 14 — Figure A.5 : Effect of the boundedness constraint ( $0\leq y\leq 1$ ) in (a) noise-free case and in (b) the presence of Gaussian process noise.

Explore similar work

PoreML: A Data-Driven Framework for Learning Multiphase Flow in Porous Media
Oct 7, 2026Chunyang Wang, Mingrui Zhang, Yuyan Zhang +5Neural Surrogate Modeling Computational Fluid Dynamics
Imperial College London · Stanford University · EarthFlow AI, Inc. +1
Discovery of Interpretable Physical Laws in Materials via Language-Model-Guided Symbolic Regression
Feb 26, 2026Yifeng Guan, Chuyi Liu, Dongzhan Zhou +4Scientific Discovery Symbolic Regression
Shanghai Artificial Intelligence Laboratory, Shanghai, 200232, China. · School of Physics, Zhejiang University, Hangzhou, 310058, China. · Soochow Institute for Energy and Materials Innovations (SIEMIS), Soochow University, Suzhou, 215006, China. +3
Physics-guided machine learning for saturation control across variable hydraulic conditions in heap leaching
Oct 11, 2026Sultan Alhamdan, Daniel Navia, Paulina QuintanillaPDE Surrogate Modeling Physics-Informed ML
Department of Chemical Engineering, University College London, London, United Kingdom · Department of Chemical and Environmental Engineering, Universidad Técnica Federico Santa María, Valparaíso, Chile

Oct 7, 2026J/K move · Enter open · S save

PoreML: A Data-Driven Framework for Learning Multiphase Flow in Porous Media

Chunyang Wang, Mingrui Zhang, Yuyan Zhang, Linqi Zhu+4

Imperial College London · Stanford University · EarthFlow AI, Inc.+1

Neural Surrogate Modeling Computational Fluid Dynamics

Open paper AskPDF arXiv

Multiphase flow in porous microstructures is central to CO

_2

storage, fuel-cell operation, and flip-chip packaging. Predicting these flows remains challenging because wettability and complex pore geometry govern the nonlinear evolution of fluid interfaces. Machine learning holds substantial promise for advancing the field, but progress is constrained by scarce time-resolved 3D datasets and a lack of a unified workflow for training and evaluating models. To fill this critical gap, we introduce PoreML, an open-source framework unifying data generation, model training, and evaluation grounded in pore-scale physics. The framework comprises three core components. (a) A modern GPU-native lattice Boltzmann solver, validated against analytical solutions and published experiments, enables reproducible data generation. (b) A 3.3 TB dataset contains 560 simulation runs and 158,546 stored time steps across four application-driven scenarios. These trajectories span synthetic structures and geometries derived from micro-CT scans of real materials, covering diverse wetting conditions and viscosity ratios. (c) A unified learning framework evaluates one-step prediction and autoregressive rollouts. Its domain-specific evaluation protocols assess predictive accuracy and physical consistency. We evaluate five models of diverse architecture under these protocols. Two complementary challenges assess transfer to larger domains and from synthetic to micro-CT-derived structures. PoreML provides a shared foundation for machine-learning research on multiphase flow in porous media, with the aim of empowering the community to develop reliable predictive models and advance the field.

Abstract

Explore similar work

PoreML: A Data-Driven Framework for Learning Multiphase Flow in Porous Media

Discovery of Interpretable Physical Laws in Materials via Language-Model-Guided Symbolic Regression

Physics-guided machine learning for saturation control across variable hydraulic conditions in heap leaching