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,Sw) space and the red axes indicating the mapped (s∗,Sw∗) space. The upper left corner of the red box corresponds to the reference point (smin,Sw,max) and the lower right corner to (sres,Sw,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∗∈[0,1] and Sw∗∈[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.
Figure 3: Illustration of (a) mutation and (b) crossover operations in an evolutionary symbolic regression algorithm.
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.
Figure 5: Data and physics loss components and the number of modes exhibited by PCSR with Lmode at varying complexity levels that optimally fit the data presented in Figure 4 .
Figure 6: Learning curves of loss components and predicted number of modes for the PCSR without and with Lmode for the data presented in Figure 4 .
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.
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.
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 .
Figure A.1 : Effect of the monotonicity constraint ( dy/dx≥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.
Figure A.2 : Effect of the monotonicity constraint ( dy/dx≥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.
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.
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.
Figure A.5 : Effect of the boundedness constraint ( 0≤y≤1 ) in (a) noise-free case and in (b) the presence of Gaussian process noise.
Multiphase flow in porous microstructures is central to CO2 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.
Chunyang Wang, Mingrui Zhang, Yuyan Zhang +5
Imperial College London · Stanford University · EarthFlow AI, Inc. +1
Discovering interpretable physical laws from high-dimensional data is a fundamental challenge in scientific research. Traditional methods, such as symbolic regression, often produce complex, unphysical formulas when searching a vast space of possible forms. We introduce a framework that guides the search process by leveraging the embedded scientific knowledge of large language models, enabling efficient identification of physical laws in the data. We validate our approach by modeling key properties of perovskite materials. Our method mitigates the combinatorial explosion commonly encountered in traditional symbolic regression, reducing the effective search space by a factor of approximately 105. A set of novel formulas for bulk modulus, band gap, and oxygen evolution reaction activity are identified, which not only provide meaningful physical insights but also outperform previous formulas in accuracy and simplicity.
Yifeng Guan, Chuyi Liu, Dongzhan Zhou +4
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
Heap-leaching hydraulics vary between regions of an ore bed, motivating fast prediction models that accommodate different material properties. This study develops a parameter-conditioned physics-guided machine-learning (PGML) surrogate and demonstrates its use in simulated saturation control. A dataset of 4096 numerical Richards-equation trajectories was divided into 2867 training, 614 validation and 615 evaluation cases. The network predicts hourly saturation profiles from the current profile, irrigation rate and four hydraulic parameters, and is trained with an approximate Gardner mass-balance penalty. A daily predictive controller compared 200 constant irrigation rates over a 24-hour horizon. Across the evaluation cases, the median full-cycle tracking deviation was 4.9% of a 0.70 saturation target, equal to that of a tuned proportional-integral (PI) controller, while the PGML controller reached the target without overshoot in the median and applied approximately 4% less irrigation. Moderate physics guidance reduced long-horizon prediction error by approximately 30% relative to data-only training and reduced irrigation-bound activity, but slightly worsened closed-loop tracking. For a matched prediction workload, the surrogate was approximately 70 times faster than the numerical solver, and batched evaluation cost under 2 ms per column for 5000 columns. The results support a reusable physics-guided surrogate for large-scale saturation control in simulation.
Sultan Alhamdan, Daniel Navia, Paulina Quintanilla
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