cs.LGOct 7, 2026
SavePoreML: A Data-Driven Framework for Learning Multiphase Flow in Porous Media
Organizations: Imperial College London · Stanford University · EarthFlow AI, Inc. · Queen Mary University of London
Abstract
Multiphase flow in porous microstructures is central to CO 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.
Figures & tables
Figure 1: PoreML simulation campaigns. (A) Primary drainage and (B) capillary trapping for CO 2 sequestration and geological energy storage, (C) GDL water intrusion for fuel-cell design, and (D) flip-chip underfill for chip packaging. Only the fluid phase of interest is shown for clarity; translucent gray meshes depict the solid.
Figure 2: PoreML framework. (A) A 3D LBM solver turns generated and CT geometries under controlled flow conditions into time-resolved phase, pressure, and velocity fields; models predict the next frame and are scored with field errors and geometry-aware metrics. (B) Scale-Up : train on small domains, test on larger ones at the same resolution. (C) Geo-Shift : train on generated geometries, test on CT scans.
Figure 3: Representative pore-scale dynamics. Localized events are highlighted in yellow; these mechanisms are not exclusive to the displayed campaign, and events such as snap-off can also occur during drainage. (A) Primary drainage : ganglion coalescence, where interfacial tension retracts menisci and connects neighboring non-wetting bodies. (B) Capillary trapping : snap-off, where wetting films invade a constricted throat and pinch a non-wetting thread into an isolated ganglion. (C) GDL intrusion : Haines jump, an abrupt meniscus depinning and pore-filling event in a fibrous throat network. (D) Underfill : microbubble or void entrapment, where advancing encapsulant fronts around solder balls leave trapped air pockets.
Figure 4: Geometry families used by PoreML. (A) Rock-like media , each shown as one large sample ( voxels) beside two small ones ( ): (A.1) sandstone from a CT scan and generated (A.2) sphere packs, (A.3) blob media, and (A.4) polydisperse spheres, with (A.5) pore and throat radii (maximal inscribed spheres). (B) Gas-diffusion layers at and voxels: (B.1) CT scan, (B.2) generated fiber mat, and (B.3) pore and throat radii. (C) Underfill : flip-chip gaps set by solder-ball diameter, gap height, and pitch. More samples: Appendix G .
| Dataset name | Cases | Domain size (voxels) | Contact angle ( ∘ ) | Viscosity ratio | Domain geometry | Simulation protocol |
| drainage_128 | 80 generated; 80 CT | Generated blobs, polydisperse spheres, and sphere packs; CT sandstones | Piston-driven non-wetting injection at ; water-wet porous outlet plate | |||
| drainage_256 | 8 generated; 8 CT | |||||
| trapping_128 | 80 generated; 80 CT | Generated and CT rock-like media | Piston-off equilibration, then wetting-water flood at ; open, pressure-anchored outlet | |||
| trapping_256 | 8 generated; 8 CT | |||||
| gdl_128 | 80 generated; 80 CT | Generated fiber mats; CT carbon paper at 5%, 20%, and 40% PTFE loading | Piston-driven through-plane water injection at ; open outlet; no plate | |||
| gdl_256 | 8 generated; 8 CT |
Table 1: PoreML dataset summary and key parameter description. for drainage, trapping, and GDL; underfill instead follows the application convention . Contact angle is measured through the non-wetting phase for the porous-media campaigns and through the wetting encapsulant for underfill.
Figure 5: In-distribution rollout curve and transfer capability. (a) In-distribution Phase MAE along the 64-step in-distribution rollout. (b) Geo-Shift and (c) Scale-Up: step-64 phase MAE relative to the same model’s in-distribution MAE (generated-only for Geo-Shift). A ratio of 1.5 indicates 50% higher error; the dashed line at 1 indicates no change.
Appendix figures & tables35 assets
Supplementary material from the paper’s appendix.
Appendix
Figure 6: The D3Q19 lattice stencil of Eq. 8 : a rest velocity, six axial neighbors (solid), and twelve face-diagonal neighbors (dashed). Arrows are numbered in the solver’s direction ordering, with opposite directions in adjacent pairs.
| Quantity | Value |
|---|---|
| Lattice | D3Q19 (rest + 6 axial + 12 diagonal) |
| Weights | (rest / axial / diagonal) |
| Sound speed | |
| Recoloring strength | (Eq. 15 ) |
| MRT non-shear rates | (drainage, trapping, GDL); fixed (underfill) |
| Wall BC | half-way bounce-back |
Table 4: Lattice constants and default numerical parameters.
| Control variable | Definition |
|---|---|
| Capillary number | |
| Viscosity ratio | |
| Contact angle | angle of the non-wetting red phase at the wall |
| Surface tension | Laplace coefficient, |
Table 5: Mapping from lattice parameters to physical control variables.
| Case | Reference |
|---|---|
| Laplace law (§ E.1 ) | |
| Layered Poiseuille (§ E.2 ) | two-layer analytic profile |
| Contact angle, flat (§ E.3 ) | imposed and Akai et al. (2018) |
| Contact angle, curved (§ E.4 ) | sphere-on-sphere geometry |
| Washburn imbibition (§ E.5 ) | two-fluid Washburn law |
| Capillary pressure (§ E.6 ) | curvature-correct entry pressures |
Table 6: Validation suite at a glance. All cases are 3D (D3Q19).
Figure 7: Laplace-law validation. (A) A deliberately non-spherical initial slab relaxes to a sphere under the interfacial force alone (sphericity ). (B) Measured pressure jump across equilibrated 3D droplets versus inverse radius , swept over input surface tension – : one fitted line per , each linear with .
Figure 8: Layered Poiseuille validation in 3D: steady velocity profiles for viscosity ratios – (LBM dots) against the two-fluid analytic solution (lines). The shaded band marks the central low-viscosity layer; maximum relative error is at .
Figure 9: Contact-angle validation on a flat wall. Left: equilibrated sessile droplets for imposed angles – , with the measured angle annotated per panel. Right: measured versus imposed angle across the sweep. The realized angle is within of the imposed value everywhere except the extreme ( ).
Figure 10: Contact-angle validation on a curved wall for imposed (A, top row) and (B, bottom row). The first three panels of each row are snapshots of the droplet (red) relaxing on the spherical grain (gray) between the volume-matched half-shell collar and the equilibrium cap; the rightmost panel is the final midplane slice, with the dotted circle marking the analytical sphere-on-sphere solution.
Figure 11: Washburn imbibition validation. (A) Penetration length from the simulation (dots) against the two-fluid Washburn solution (lines) for contact angles – . (B) The tube at five successive times, rendered upright (wetting phase red, rising bottom to top under capillary suction alone, ): 3D view above, midplane slice below.
Figure 12: Primary-drainage capillary-pressure curves ( versus wetting saturation ) for the circular tube (top left) and the rectangular tube (top right), swept over contact angle. Dashed lines mark each angle’s curvature-correct analytic entry pressure (circle: Mason–Morrow shape-factor form on the voxelized cross-section; rectangle: two-curvature Laplace). Both tubes drain along a flat plateau that tracks the analytic entry per angle. Bottom: the non-wetting phase (red) crossing the throat and filling the duct piston-like (upper sequence: circular tube; lower sequence: rectangular tube; geometry at right).
Figure 13: Dynamic primary drainage in a four-throat micromodel. (a) Experimental images adapted from Figure 4(a) of Zacharoudiou et al. (2017) . (b) Midplane slices of our three-dimensional color-gradient LBM simulation at corresponding meniscus positions ( , imposed ). From left to right, the non-wetting phase advances through the inlet, fills the pore body, and invades the widest outlet throat, ( in the experiment, in the simulation); and remain predominantly wetting-filled. In the simulation, red denotes the non-wetting phase, white the wetting phase, and gray the solid. The dashed arrow indicates the inlet flow direction.
Figure 14: Spontaneous imbibition in a four-throat micromodel. (a) Experimental images adapted from Figure 8(a) of Zacharoudiou et al. (2017) . (b) Midplane slices of our three-dimensional color-gradient LBM simulation at corresponding meniscus positions (imposed , ). From left to right, the wetting phase enters through , advances into the pore body, and invades the adjacent throat in stage (3), filling it by stage (4), before the narrowest throat . The stage-(3) heading describes the experiment; the simulation at that stage shows partial invasion of . Their experimental widths are and , represented by and in the simulation. The remaining non-wetting phase then withdraws through the widest throat, . In the simulation, red denotes the non-wetting phase, white the wetting phase, and gray the solid. The dashed arrow indicates the inlet flow direction.
Figure 15: Displacement regimes on the sphere pack, varying only ; the non-wetting phase is red, solid translucent, flow left to right. The central map sketches the Lenormand et al. (1988) phase diagram in , and each regime region is linked to a snapshot of the corresponding run: stable displacement at — a compact front sweeping the pore space; viscous fingering at — thin branched fingers that screen the surrounding pore space; capillary fingering at — a rough, capillarity-controlled front.
| Campaign | Domain (voxels) | Runs | Frames per run | Frame interval | Frames |
|---|---|---|---|---|---|
| Drainage | 160 | 376 | 5,000 | 58,830 | |
| 16 | 694 | 5,000 | 10,643 | ||
| Trapping | 160 | 105 | 5,000 | 19,424 | |
| 16 | 239 | 5,000 | 3,616 | ||
| GDL | 160 | 292 | 5,000 | 45,437 | |
| 16 | 486 | 5,000 | 8,217 |
Table 7: Run lengths per campaign and domain size. Domains include the reservoir, plate, and buffer slabs; frames include the initial state. Frames per run are medians.
Figure 16: Input geometry examples at the large domain size. Each panel renders the solid phase of four samples of one family. (a)–(c) Sandstones from CT scans, voxels: Bentheimer, Buff Berea, and CastleGate. (d)–(f) Generated rock-like media, voxels: sphere packs, blob media, and polydisperse overlapping spheres. (g)–(i) Carbon-paper GDLs from CT scans at 5%, 20%, and 40% PTFE loading, voxels, with the through-plane flow direction vertical. Samples of one family share a morphology but differ in porosity and connectivity.
Figure 17: Specks corrupt the Euler number. (a) counts pieces, tunnels, and enclosed bubbles of the non-wetting phase, so every island of non-wetting phase and every bubble of wetting phase adds , however small: the blob alone has , with the specks . (b, c) One -voxel crop of a slice across the flow, fiber GDL, rollout step 64: the reference and the thresholded prediction ( ) of P3D and of Transolver++, with speck islands (black) and speck bubbles (yellow) marked. Under each panel, of the whole 3D frame with its three terms. The predictions’ 84 and 484 pieces and 95 and 8 bubbles are almost all specks: P3D’s has the wrong sign against the reference’s , and Transolver++’s hides 587 tunnels behind 484 islands.
Figure 18: Tunnels corrupt the Euler number and cannot be filtered. (a) A loop of the non-wetting phase is a tunnel through it and lowers by one; being part of the phase, no size filter removes it. Where the reference keeps a wetting film between a fiber and the phase there is no loop ( ); a prediction that engulfs the fiber, or joins two regions by two voxel-thin bridges, closes one ( ). (b) Reference and FNO3D in a plane along the flow, where the fibers are cut across ( voxels): FNO3D encircles three fiber cross-sections that in the reference touch wetting fluid. (c) Reference and Transolver++ in a plane across the flow: twelve pockets of wetting fluid enclosed by thin strands of phase. A region the phase encircles in the plane, a fiber cross-section or a pocket of wetting fluid, is the hole of a tunnel and is outlined in black. Under each panel, and the tunnel count of the whole frame, and the loop count of the crop.
Figure 19: Drainage, in-distribution: three-dimensional rollout. A generated blob medium at voxels with porosity , and ; the rollout starts at stored frame 128. Rows are the reference simulation and the five baselines, columns the six rollout steps, and each predicted panel is labeled with its IoU. Over the window the reference inflates a single ganglion in the center of the domain, from a bud at step 7 to a rounded body at step 64, and the non-wetting saturation rises by half, from to . AB-UPT reproduces that growth at IoU and ends on the reference saturation; UNet3D ( ) and P3D ( ) follow it slightly full; FNO3D ( ) over-invades to , and Transolver++ ( ) fragments the interface as the rollout proceeds.
Figure 20: Drainage, in-distribution: mid-plane slice. The same run and steps as Figure 19 , sliced through the middle of the domain. The inflating ganglion lies off this plane, so the reference slice holds its shape for the whole window — which makes the figure a test of what a rollout invents where nothing happens. Every model eventually opens a spurious non-wetting patch inside the wetting-filled pore below the main body: Transolver++ carries one from step 18 to the end, FNO3D grows a large one from step 53, and UNet3D and P3D each open one at step 64.
Figure 21: Capillary trapping, in-distribution: three-dimensional rollout. A generated blob medium at voxels with porosity , and ; the rollout starts at stored frame 64. The wetting flood removes nearly a third of the trapped phase, from to saturation, and breaks the connected body in the middle of the domain into separate ganglia. Every baseline under-mobilizes: the endpoint saturations run from (UNet3D) to (FNO3D) against the reference , and the step-64 overlaps spread from (UNet3D), (FNO3D) and (AB-UPT) down to (Transolver++) and (P3D), which keeps the central body nearly intact.
Figure 22: Capillary trapping, in-distribution: mid-plane slice. The same run and steps as Figure 21 ; here red is the invading wetting phase and blue the non-wetting phase being displaced. The reference thins the trapped blob at the top of the slice from the side facing the invading phase and narrows the neck that connects it to the throat on its left. All five models keep the blob fuller than the reference; P3D keeps it widest and its overlap decays furthest, from at step 7 to at step 64.
Figure 23: GDL intrusion, in-distribution: three-dimensional rollout. A generated fiber mat at voxels with porosity , and ; the rollout starts at stored frame 64. Water climbs through the mat from the inlet at the bottom and the saturation doubles, from to — the largest relative change of any run in this appendix. The baselines bracket that growth from both sides: UNet3D ends half again too wet at , FNO3D at , AB-UPT closest at , while P3D stops at and Transolver++ at , barely above the saturation it started from. The step-64 overlaps follow the same order: FNO3D , UNet3D , AB-UPT , P3D , Transolver++ .
Figure 24: GDL intrusion, in-distribution: mid-plane slice. The same run and steps as Figure 23 . In the slice the invading phase goes from two small patches at the bottom to a single body spanning the middle of the mat. The errors are errors of timing: UNet3D, P3D and AB-UPT merge the two patches into one cluster already at step 7, ahead of the reference, while Transolver++ starts correctly and then stops growing, leaving a thin ragged filament that never fills the pore space — the slice-level counterpart of its endpoint saturation.
Figure 25: Drainage, Geo-Shift: three-dimensional rollout. A CastleGate sandstone sample at voxels with porosity , and ; the rollout starts at stored frame 128. The reference gains almost half its non-wetting volume, from to saturation, extending the cluster along the pore and opening a new branch late in the window. How much of that branch a model produces separates the rollouts: AB-UPT reaches IoU and is the only baseline to end below the reference saturation, at ; UNet3D and P3D end just above it; Transolver++ roughens the interface, and FNO3D over-invades to .
Figure 26: Drainage, Geo-Shift: mid-plane slice. The same run and steps as Figure 25 . The slice isolates the new branch: the arm at the top of the slice is wetting-filled at step 7 and half non-wetting by step 64. UNet3D, P3D and AB-UPT advance that front at about the right rate, Transolver++ carries it higher up the arm than the reference, and FNO3D additionally opens a pore on the right that the reference leaves filled with the wetting phase — which is where its saturation excess accumulates.
Figure 27: Capillary trapping, Geo-Shift: three-dimensional rollout. A Bentheimer sandstone sample at voxels with porosity , and ; the rollout starts at stored frame 64. The wetting flood strips about a fifth of the trapped non-wetting phase, from to saturation, and the models bracket that endpoint: P3D and UNet3D remove too much and end at , FNO3D and Transolver++ too little at and , and AB-UPT lands closest at with the highest overlap, .
Figure 28: Capillary trapping, Geo-Shift: mid-plane slice. The same run and steps as Figure 27 ; red is the invading wetting phase and blue the non-wetting phase being displaced. The wetting phase eats into the upper body from the right, and where a model stops that erosion sets its endpoint. The overlap need not decay monotonically: UNet3D dips to at step 41, when it has cleared too much of the body, and recovers to by step 64 as the reference catches up with it.
Figure 29: GDL intrusion, Geo-Shift: three-dimensional rollout. A CT carbon-paper GDL with 40% PTFE coating at voxels, porosity , and ; the rollout starts at stored frame 64. Water invasion doubles over the window, from to saturation, and the baselines bracket that endpoint widely: UNet3D and FNO3D end half again too wet, at and , P3D lands on the reference volume at , and Transolver++ stops short at . The step-64 overlaps are P3D , AB-UPT , UNet3D and FNO3D , Transolver++ .
Figure 30: GDL intrusion, Geo-Shift: mid-plane slice. The same run and steps as Figure 29 . In the scanned layer the water advances as a broad body between the fibers rather than as a front, and the models disagree about where it starts: UNet3D and P3D nucleate a large cluster at the lower left already at step 7, FNO3D builds one at the right edge that the reference reaches only at step 53, and Transolver++ leaves isolated specks that never merge, its overlap decaying to while its saturation stalls at .
Figure 31: Drainage, Scale-Up: three-dimensional rollout. A Buff Berea sandstone sample at voxels with porosity , and ; the rollout starts at stored frame 192. In eight times the training volume the reference extends its cluster through the pore network and moves from to saturation. FNO3D holds the pattern best here at IoU , just ahead of AB-UPT , which also ends closest in volume; P3D and UNet3D keep the shape but end too wet, at and ; Transolver++ shreds the body into disconnected fragments and ends at saturation, half again the reference.
Figure 32: Drainage, Scale-Up: mid-plane slice. The same run and steps as Figure 31 . The slice shows where the excess volume accumulates: UNet3D and P3D thicken the cluster where it meets the wetting phase at the top, and Transolver++ breaks into the wetting-filled channel at the lower left from step 30 onward and holds that intrusion to the end of the window, while the reference keeps the channel filled.
Figure 33: Capillary trapping, Scale-Up: three-dimensional rollout. A generated sphere pack at voxels with porosity , and ; the rollout starts at stored frame 64. The flood removes a sixth of the trapped phase, from to saturation, while the many isolated ganglia stay where they are: at this size the window covers a slow, distributed erosion rather than a single event. Overlap and volume disagree here. Transolver++ holds the highest overlap, , while ending furthest above the reference saturation at , and UNet3D ends within of the reference volume at the lowest overlap, ; AB-UPT alone strips too much and ends at .
Figure 34: Capillary trapping, Scale-Up: mid-plane slice. The same run and steps as Figure 33 ; red is the invading wetting phase and blue the non-wetting phase being displaced. The slice makes the slowness explicit: the ganglion at the center keeps its outline for the whole window and loses volume only along its rim, and the models differ by a few voxels around that rim — AB-UPT erodes it furthest, Transolver++ least.
Figure 35: GDL intrusion, Scale-Up: three-dimensional rollout. A generated fiber mat at voxels with porosity , and ; the rollout starts at stored frame 192. The reference invades a third more of the layer, from to saturation, across four times the in-plane area of the training domain. P3D ( ) and UNet3D ( ) hold the intrusion pattern best; FNO3D and Transolver++ (both ) carry too much water and end at against the reference ; AB-UPT ends on the right volume, , at the lowest overlap of the five, .
Figure 36: GDL intrusion, Scale-Up: mid-plane slice. The same run and steps as Figure 35 . In the slice the water holds a flat band across the bottom of the layer and thickens it instead of fingering upward, and the models differ along the right edge of that band: AB-UPT extends it across the foot of the slice from step 30 onward, UNet3D and FNO3D follow later, and P3D keeps the band closest to the reference width.