Physics-informed convolutional neural networks for fluid flow through porous media

# Physics-informed convolutional neural networks for fluid flow through porous media
Source: [https://arxiv.org/html/2605.20250](https://arxiv.org/html/2605.20250)
###### Abstract

Accurate simulation of fluid flow in porous media is a challenging task due to the complexity of pore\-space geometries and the computational cost of solving the Navier–Stokes equations\. Traditional numerical solvers rely on carefully constructed meshes, often requiring manual intervention, and typically exhibit slow convergence\. This difficulty is particularly pronounced in porous media, where the diffusive nature of momentum transport is hindered by intricate solid boundaries\. These challenges limit the efficiency of numerical simulations, particularly when repeated evaluations are required\. We present a neural\-network\-based framework for predicting pore\-scale velocity fields directly from sample geometry\. The method is based on a convolutional encoder–decoder architecture with skip connections, designed to preserve fine\-grained structural information\. Physical consistency is encouraged through a custom loss function composed of multiple terms: incompressibility, no\-flow conditions within solids, periodicity constraints, and agreement with the global tortuosity index\. We systematically analyze the influence of weight selection for these loss terms, quantifying their individual contributions to prediction accuracy\. Several architectural variants inspired by computer vision are evaluated to identify one providing the best performance and robustness\. The generalization ability of the trained network is assessed on samples outside the training distribution, including variations in boundary conditions, obstacle geometry, and porosity\. Finally, we demonstrate additional practical applications in which network predictions are used to initialize the Lattice–Boltzmann simulations, a standard fluid dynamics solver, often used in complex boundary problems like porous media and used by us to train the network\. We have used network\-generated velocity field as a starting point and found that this significantly accelerates LBM solver convergence, achieving improvements in over 90% of cases\.

###### keywords:

physics\-informed neural networks , pore\-scale velocity prediction , convolutional encoder–decoder networks , model generalization , fluid flow , porous media , tortuosity and permeability , the Lattice Boltzmann Method

††journal:Scientific Reports\\affiliation

organization=Dioscuri Center in Topological Data Analysis, Institute of Mathematics, Polish Academy of Sciences,addressline=ul\. Śniadeckich 8,postcode=00\-656,city=Warsaw,country=Poland\\affiliationorganization=Institute of Experimental Physics, Faculty of Physics and Astronomy, University of Wrocław, addressline=pl\. M\. Borna 9,postcode=50\-204,city=Wrocław,country=Poland

\\affiliation

organization=Dioscuri Center in Topological Data Analysis, Institute of Mathematics, Polish Academy of Sciences,addressline=ul\. Śniadeckich 8,postcode=00\-656,city=Warsaw,country=Poland

\\affiliation

organization=Institute of Theoretical Physics, Faculty of Physics and Astronomy, University of Wrocław,addressline=pl\. M\. Borna 9,city=Wrocław,postcode=50\-204,country=Poland\\affiliationorganization=Parallel and Distributed Systems Laboratory, Jožef Stefan Institute,addressline=Jamova cesta 39,city=Ljubljana,postcode=1000,country=Slovenia\\affiliationorganization=\*Corresponding author: maciej\.matyka@uwr\.edu\.pl \(Maciej Matyka\)

## 1Introduction

Porous materials are ubiquitous in nature and technology\. Research on them covers a wide range of topics related to fluid transport, mechanics, and strength\. One area with particular applications is the study of the permeability of porous samples, which has applications in oil extractionZhong et al\. \([2021](https://arxiv.org/html/2605.20250#bib.bib37)\), CO2storageCossins et al\. \([2023](https://arxiv.org/html/2605.20250#bib.bib3)\), and in medicine, with a prominent example of the blood\-brain barrier as the drug delivery path, of which permeability is often of interestFong \([2015](https://arxiv.org/html/2605.20250#bib.bib6)\)\.

For the past few decades, solving the fluid transport problem required the use of advanced numerical methods to solve the Navier\-Stokes equationsAnderson \([2002](https://arxiv.org/html/2605.20250#bib.bib2)\)\. The main time\-demanding elements of standard computational fluid dynamics simulation pipelines are: solving large systems of linear equations, creating computational grids, utilizing high\-powered computers, and the diffusive nature of momentum transport\. In general, given boundary conditions, initial conditions, and physical parameters, the only way to find the trajectories of particles in water is to solve the non\-linear Navier\-Stokes equations in approximate, numerical, and iterative approaches\. The use of a mesoscopic Lattice\-Boltzmann MethodGuo and Shu \([2013](https://arxiv.org/html/2605.20250#bib.bib12)\), in which we deal with the transport of a particle distribution function on regular grids, is promising and successful, particularly in irregular and porous geometriesSucci \([2001](https://arxiv.org/html/2605.20250#bib.bib32)\)\. It comes, however, at an additional physical complexity in interpretation and additional cost of memory and time demand, thus requiring special optimizationsLehmann et al\. \([2022](https://arxiv.org/html/2605.20250#bib.bib22)\)and memory layouts for handling large datasetsTomczak and Szafran \([2019](https://arxiv.org/html/2605.20250#bib.bib34)\)\.

Meanwhile, we \(humans\) learn by observing physical systems\. Looking at the ball thrown to the sky, we know that it will eventually fall, and we do this naturally, not knowing the PDEs behind\. Similarly, for fluid flow in an open channel with an obstacle, we know approximately where the water will go and turn, as well as under what conditions a vortex may form\. It is intuitively consistent that vortex\-dominated flow occurs under conditions of increased velocity, and we are able to predict based on our knowledge and previous observations\. Recent research has shown that this is also possible using artificial intelligence and deep convolutional neural networks\. Guo et al\. have demonstrated steady\-state flow predictions for the flow over single objects immersed in a channel, utilizing convolutional neural networks \(CNNs\)Guo et al\. \([2016](https://arxiv.org/html/2605.20250#bib.bib11)\)\. Recently, the same authors have shown unsteady flow predictions over single obstacles in a channelGuo et al\. \([2024](https://arxiv.org/html/2605.20250#bib.bib10)\)\. We are, however, aiming to predict the flow in complex geometries of porous media, both at low and high porosity, where the complexity of the pore network poses a challenge for neural networks\. In our recent work we have shown, that only geometry is enough to predict physical properties of porous samples to predict their macroscopic properties in fluid flow and diffusive processes using CNN’sGraczyk and Matyka \([2020](https://arxiv.org/html/2605.20250#bib.bib8)\); Graczyk et al\. \([2023](https://arxiv.org/html/2605.20250#bib.bib9)\)\. Research record on the topic of predicting physical properties of physical samples using deep learning techniques is rich\. Recently, Lin et al\. studied the application of CNNs in predicting velocity and temperature distributions in metal foam samples, demonstrating the outstanding efficiency of the method compared to standard numerical proceduresLin et al\. \([2025](https://arxiv.org/html/2605.20250#bib.bib24)\)\. It has been shown that the prediction of physical fields using CNNs can be improved by incorporating a physically informed network to reproduce temperature fieldsZhao et al\. \([2023](https://arxiv.org/html/2605.20250#bib.bib36)\), a procedure that can be equivalent to solving differential equations and serves to generate physically correct results\.

Our work proposes using a combined physics\-informed CNN approach for fluid flow in complex porous media at various porosity conditions\. We use convolutional neural networks and a custom, physically motivated loss function\. The main contributions of our work within CNN\-based pore\-scale flow prediction are summarized as follows:

- 1\)We implemented a complete process from constructing the porous samples, simulation using the Lattice Boltzmann Method, through CNN training for predicting velocity fields in fluid flow through a porous media at various porosities\.
- 2\)We construct and evaluate a custom loss function tailored to pore\-scale flow prediction, including terms related to periodicity of the computational grid and tortuosity matching\. \.
- 3\)We carried out tests on various types of porous media samples to demonstrate the generalization abilities of the neural network and predict beyond the data region on which the network was learned\.
- 4\)We show practical application of CNNs for the Lattice Boltzmann Method: LBM solver can be significantly accelerated if the neural network’s predictions are used to initialize the fluid dynamics solver\.

We note that the general use of CNNs as surrogate models and the idea of incorporating physics\-inspired constraints into neural\-network training are well established\. The contribution of the present work is therefore not the introduction of a new learning paradigm, but rather the formulation and evaluation of a problem\-specific physics\-informed loss for pore\-scale flow prediction, together with a systematic assessment of its components, architecture dependence, out\-of\-distribution behavior, and use for LBM initialization\.

The paper is organized as follows: in Sec\.[2](https://arxiv.org/html/2605.20250#S2), we introduce the physical model of porous media and discuss the fluid flow solver, the Lattice Boltzmann Method, which is used later to prepare the training velocity fields\. We describe the procedure to synthesize porous medium models investigated in the paper\. Additionally, the learning process is described in detail, including data augmentation and a detailed description of the loss function\. We discuss the choice and analysis of the neural network architecture and describe details of the training protocol\. In Sec\.[Results](https://arxiv.org/html/2605.20250#Sx1)we present the results of predictions made by the trained network, both for tortuosity and permeability of the samples\. Then, we discuss and present results on the generalization properties of our approach, including its ability to generalize to other porosities, shapes of obstacles, and other boundary conditions\. We summarize our findings in Sec\.[3](https://arxiv.org/html/2605.20250#S3)\.

## 2Materials and Methods

### Pore scale flow

We describe porous samples at the pore\-scale level, where solid obstacles, distributed randomly, fill the space up to the desired porosity of the sample\. We assume periodic boundary conditions in both directions, which, for the two\-dimensional case, directly represent a porous medium on a torus topology\. The flow is driven by gravity, which accelerates the fluid to a certain level\. Due to the no\-slip boundary conditions, the system stabilizes to a so\-called stationary flow\. The governing equations for the fluid flow in pores are Navier\-Stokes, that, for an incompressible fluid takes form:

∇⋅𝐮=0,\\nabla\\cdot\\mathbf\{u\}=0,\(1\)ϱ​\(𝐮⋅∇𝐮\)=−∇p\+μ​∇2𝐮\+𝐟,\\varrho\(\\mathbf\{u\}\\cdot\\nabla\\mathbf\{u\}\)=\-\\nabla p\+\\mu\\nabla^\{2\}\\mathbf\{u\}\+\\mathbf\{f\},\(2\)where𝐮\\mathbf\{u\}is velocity,μ\\muis dynamic viscosity,ppis the pressure, andϱ\\varrhois its density\. Here, we operate in the regime of low Reynolds numbers, for which the flow through porous matrix satisfies simpler Darcy’s linear law:

𝐪=−kμ​\(∇P\+ϱ​𝐟\),\\mathbf\{q\}=\-\\frac\{k\}\{\\mu\}\\left\(\\nabla P\+\\varrho\\mathbf\{f\}\\right\),\(3\)whereqqis Darcy flux,kkis permeability\[L2\]\[L^\{2\}\]andffis the external force density \(i\.e\. gravity\)\. Permeability,kk, is of utmost importance in applications and porous media research\. It became a standard measure with a described protocol for the experimental treatment\. For practical applications, it can be expressed in terms of tortuosity and porosity, i\.e\.:

k=φ3c​τ2​S,k=\\frac\{\\varphi^\{3\}\}\{c\\tau^\{2\}S\},\(4\)whereccis a material\-specific constant,τ\\tauis tortuosity of the medium, andSSis specific surface area of poresKoponen et al\. \([1997](https://arxiv.org/html/2605.20250#bib.bib19)\)\. Such a definition of permeability allows for estimations based on the geometrical properties of the medium, except for the value ofτ\\tau, which must be measured in the simulation or experiment\. Tortuosity is a non\-dimensional number representing elongation of porous media pores due to existence of solid matrix and it is driven by pore space geometry and physical process of channel formation\. It can be investigated using streamlines generated on top of the velocity field, whose averaged elongation can be used to compute the tortuosity indexKoponen et al\. \([1996](https://arxiv.org/html/2605.20250#bib.bib18)\); Matyka et al\. \([2008](https://arxiv.org/html/2605.20250#bib.bib26)\)\. We will use the velocity field to compute tortuosity by directly integrating the streamwise and total momentum in the pore space, which results in the following expression for tortuosity:

τ=⟨v⟩⟨vx⟩,\\tau=\\frac\{\\langle v\\rangle\}\{\\langle v\_\{x\}\\rangle\},\(5\)where brackets denote average over pore space, andxxis the component of velocity parallel to the direction of external body forceDuda et al\. \([2011](https://arxiv.org/html/2605.20250#bib.bib5)\); Matyka and Koza \([2012](https://arxiv.org/html/2605.20250#bib.bib27)\)\.

### The Lattice Boltzmann Method

To solve the flow problem and obtain pore\-space velocity fields to train the neural network, we utilize an indirect, mesoscopic Lattice Boltzmann Method \(LBM\)\. LBM works on a regular grid with a nine\-component velocity discretization\. It solves the Navier\-Stokes flow equations indirectly by solving the transport equation of the particle distribution functionf​\(𝐱,t\)f\(\\mathbf\{x\},t\), which represents the probability of finding particles at a given velocity at a specific location and time\. The method demonstrated its unique properties and accuracy in porous media flows, becoming popular primarily due to its simplicity in implementation and the ability to handle complex geometries using local no\-slip boundary conditions\. We will utilize its properties and solve the transport equation of the particle distribution function using the following model in the BGK \(one\-relaxation approximation\) collision termSucci \([2001](https://arxiv.org/html/2605.20250#bib.bib32)\):

fi​\(𝐱\+𝐞i,t\+δ​t\)=fi​\(𝐱,t\)\+fie​q​\(𝐱,t\)−fi​\(𝐱,t\)τ,f\_\{i\}\(\\mathbf\{x\}\+\\mathbf\{e\}\_\{i\},t\+\\delta t\)=f\_\{i\}\(\\mathbf\{x\},t\)\+\\frac\{f\_\{i\}^\{eq\}\(\\mathbf\{x\},t\)\-f\_\{i\}\(\\mathbf\{x\},t\)\}\{\\tau\},\(6\)whereτ\\tauis relaxation time, andfie​qf\_\{i\}^\{eq\}is distribution function in equilibrium\. Calculation of consecutive moments of the transported distribution function gives the macroscopic densityϱ​\(𝐱,t\)\\varrho\(\\mathbf\{x\},t\)and velocity𝐮​\(𝐱,t\)\\mathbf\{u\}\(\\mathbf\{x\},t\)fields:

ϱ​\(𝐱,t\)=∑ifi​\(𝐱,t\),\\varrho\(\\mathbf\{x\},t\)=\\sum\_\{i\}f\_\{i\}\(\\mathbf\{x\},t\),\(7\)ϱ​\(𝐱,t\)​𝐮​\(𝐱,t\)=∑ifi​\(𝐱,t\)​𝐞𝐢,\\varrho\(\\mathbf\{x\},t\)\\mathbf\{u\}\(\\mathbf\{x\},t\)=\\sum\_\{i\}f\_\{i\}\(\\mathbf\{x\},t\)\\mathbf\{e\_\{i\}\},\(8\)where𝐞𝐢\\mathbf\{e\_\{i\}\}are the lattice vectors\.

### Porous samples

Two\-dimensional porous structures were represented as binary images with a resolution of256×256256\\times 256pixels\. The structures were generated by superimposing standing sinusoidal waves with a fixed wavelength and amplitude but varying random directions𝐪i\\mathbf\{q\}\_\{i\}and phasesϕi\\phi\_\{i\}\. The resulting field is given by:

f​\(𝐱\)=2N​∑i=1Ncos⁡\(𝐪i⋅𝐱\+ϕi\),𝐱∈\[0,1\]×\[0,1\],f\(\\mathbf\{x\}\)=\\sqrt\{\\frac\{2\}\{N\}\}\\sum\_\{i=1\}^\{N\}\\cos\(\\mathbf\{q\}\_\{i\}\\cdot\\mathbf\{x\}\+\\phi\_\{i\}\),\\quad\\mathbf\{x\}\\in\[0,1\]\\times\[0,1\],\(9\)
where𝐱\\mathbf\{x\}is the position vector, and the phasesϕi\\phi\_\{i\}are randomly sampled from the interval\[0,π\]\[0,\\pi\]\. The wave vectors are defined as𝐪i,j=2​k​π\\mathbf\{q\}\_\{i,j\}=2k\\pi, wherekkis a randomly chosen integer within the range\[−12,12\]\[\-12,12\]\. For each generated structure, the number of superimposed standing waves,NN, is randomly selected beforehand from the range\[10,100\]\[10,100\]\.

Given the random function \([9](https://arxiv.org/html/2605.20250#S2.E9)\), the different phases of the system are defined as

𝐱∈ℬ\\displaystyle\\mathbf\{x\}\\in\\mathcal\{B\}iff​\(𝐱\)≤ε,\\displaystyle f\(\\mathbf\{x\}\)\\leq\\varepsilon,𝐱∈𝒫\\displaystyle\\mathbf\{x\}\\in\\mathcal\{P\}iff​\(𝐱\)\>ε\.\\displaystyle f\(\\mathbf\{x\}\)\>\\varepsilon\.\(10\)Here,ℬ\\mathcal\{B\}represents the solid phase, while𝒫\\mathcal\{P\}corresponds to the pore phase\. The threshold valueε\\varepsilonis chosen to achieve the desired porosityφ\\varphifor each generated sample\. The porosity values are uniformly sampled within the range0\.70≤φ≤0\.950\.70\\leq\\varphi\\leq 0\.95\.

In this context, the porosity is defined as the ratio of void\-phase pixels to the total number of pixels\. Structures that do not allow fluid flow through the material were removed from the dataset\. A total of 10,000 structures were generated and used for model training\.

Example of structures generated using this algorithm, along with its corresponding velocity fields, are presented in Figure[1](https://arxiv.org/html/2605.20250#S2.F1)\.

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/lbm0.png)

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/lbm1.png)

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/lbm2.png)

Figure 1:Example porous structures and corresponding LBM solutions \- color indicates velocity magnitude\. The same color scale applied on all plots\. Structure porosity is \(a\) 0\.771, \(b\) 0\.861, \(c\) 0\.917\.![Refer to caption](https://arxiv.org/html/2605.20250v1/x1.png)Figure 2:Statistical summaries of the LBM solutions for the dataset used to train the model\. \(a\) Distribution of the trivial porosity \(defined as the ratio of void area to total sample area\)\. \(b\) Distribution of the mean and maximum fluid velocity within the porous structures\. \(c\) Distribution of permeability values\. \(d\) Distribution of tortuosity values across all samples\. Red vertical lines on panels b\) and c\) indicate the prediction accuracy as indicated by the RMSE for the best model – see Table[2](https://arxiv.org/html/2605.20250#Sx1.T2)\.Figure[2](https://arxiv.org/html/2605.20250#S2.F2)presents the distributions of key physical and geometric properties derived from the dataset used for model training\.

### Data augmentations

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/aug_all.png)Figure 3:a\) Example porous structure \(top row\) and corresponding LBM solution \- color indicates velocity magnitude\. b\) The same structure and velocity field after vertical flip is applied\. c\) Effect of horizontal and vertical roll: 20% of image rows and columns \(indicated by dashed lines in panel a\) where shifted beyond the original image and placed at the other end\.During training, data augmentation techniques are applied to the structures to enhance the model’s generalization ability\. The augmentations include random vertical flipping and random rolling\. In the vertical flip augmentation, each structure is flipped along the vertical axis with a given probability \(set to0\.50\.5in all experiments\)\. The corresponding velocity field is adjusted accordingly by inverting the vertical velocity component \(vy→−vyv\_\{y\}\\rightarrow\-v\_\{y\}\), see Figure[3](https://arxiv.org/html/2605.20250#S2.F3)b\. The random roll augmentation shifts a selected number of columns from the left side of the structure to the right and similarly moves a number of top rows to the bottom\. The velocity field is adjusted consistently to maintain continuity\. The effect is shown in Figure[3](https://arxiv.org/html/2605.20250#S2.F3)c\. The number of shifted columns and rows is randomly selected for each instance, ensuring that no more than 30% of the image is rolled\.

These augmentations allow the model to learn from a more diverse set of training samples, effectively increasing the variability of the dataset without requiring additional computationally expensive LBM simulations\.

### Model architecture

The solution employs a U\-Net architectureRonneberger et al\. \([2015](https://arxiv.org/html/2605.20250#bib.bib30)\), a widely adopted neural network originally developed for pixel\-wise image segmentation\. U\-Net features an encoder\-decoder structure with skip connections, which facilitate the preservation of spatial details while enabling the extraction of hierarchical, multi\-scale features\. The input to the model is a256×256256\\times 256binary image representing the porous structure\. The image is processed through the encoder path, which is composed of a deep CNN backbone\. In our study, we consider a range of architectures commonly used in computer vision, including classical models such as VGGSimonyan and Zisserman \([2015](https://arxiv.org/html/2605.20250#bib.bib31)\), ResNetHe et al\. \([2015](https://arxiv.org/html/2605.20250#bib.bib13)\), and DenseNetHuang et al\. \([2018](https://arxiv.org/html/2605.20250#bib.bib15)\), as well as more recent designs such as ConvNeXtLiu et al\. \([2022](https://arxiv.org/html/2605.20250#bib.bib25)\), EfficientNetTan and Le \([2019](https://arxiv.org/html/2605.20250#bib.bib33)\), and MobileNetV3Howard et al\. \([2019](https://arxiv.org/html/2605.20250#bib.bib14)\)\. These newer architectures have been proposed as either modernized or computationally efficient variants of standard convolutional networks\. In particular, ConvNeXt incorporates several design elements inspired by transformer\-based models while remaining fully convolutional, EfficientNet introduces a compound scaling strategy, and MobileNetV3 focuses on lightweight computation through the use of depthwise separable convolutions and optimized building blocks\. All these architectures are used interchangeably as encoder backbones within the same framework, allowing for a consistent comparison\.

Next, a bottleneck layer refines the internal feature representation before upsampling\. The decoder consists of transposed convolutional layers that successively restore spatial resolution\. Additionally, skip connections from corresponding encoder layers are incorporated, mitigating information loss due to downsampling\. These skip connections enhance the network’s ability to produce pixel\-wise velocity predictions\.

The final output is a256×256256\\times 256velocity field with two channels, corresponding to thexxandyycomponents of the predicted velocity field \(uxu\_\{x\}anduyu\_\{y\}\)\. This prediction is compared to the reference velocity field obtained via LBM simulation, and the discrepancy is quantified using the custom loss function defined in Equation \([11](https://arxiv.org/html/2605.20250#S2.E11)\)\. The loss is then backpropagated through the network to update the model parameters\.

### Learning procedure

For each porous structure in the dataset, the LBM solution was computed\. The solution at each pixel is represented as a two\-dimensional velocity vector𝐮=\(ux,uy\)\\mathbf\{u\}=\(u\_\{x\},u\_\{y\}\)\. During the network training, the model receives the structure \(binary image\) as input and the corresponding LBM solution \(pixel\-size vector field\) as output\. The goal of training is to predict the LBM velocity field exclusively from the provided structure image of the structure\.

A notable feature of our approach is that the model backbone, which we use to construct the U\-Net architecture, consists of convolutional neural networks \(CNNs\) of different architectures\. Since CNNs operate as pattern recognition models, they are not explicitly aware of solving a differential equation\. Instead, the network learns a mapping between the input structure and the output velocity field\. To train the model, we employ a custom loss functionℒ\\mathcal\{L\}, which is composed of multiple terms:

ℒ\\displaystyle\\mathcal\{L\}=ℒvel\+α​ℒobstacle\+β​ℒdiv\+γ​ℒperio\+δ​ℒtort\\displaystyle=\\mathcal\{L\}\_\{\\text\{vel\}\}\+\\alpha\\mathcal\{L\}\_\{\\text\{obstacle\}\}\+\\beta\\mathcal\{L\}\_\{\\text\{div\}\}\+\\gamma\\mathcal\{L\}\_\{\\text\{perio\}\}\+\\delta\\mathcal\{L\}\_\{\\text\{tort\}\}\(11\)=12​\|Ω\|​∑j∈\{x,y\}∫Ω\(vj​\(𝐱\)−uj​\(𝐱\)\)2​𝑑𝐱\+α2​\|ℬ\|​∑j∈\{x,y\}∫ℬ\|uj​\(𝐱\)\|​𝑑𝐱\\displaystyle=\\frac\{1\}\{2\|\\Omega\|\}\\sum\_\{j\\in\\\{x,y\\\}\}\\int\_\{\\Omega\}\\left\(v\_\{j\}\(\\mathbf\{x\}\)\-u\_\{j\}\(\\mathbf\{x\}\)\\right\)^\{2\}d\\mathbf\{x\}\+\\frac\{\\alpha\}\{2\|\\mathcal\{B\}\|\}\\sum\_\{j\\in\\\{x,y\\\}\}\\int\_\{\\mathcal\{B\}\}\|u\_\{j\}\(\\mathbf\{x\}\)\|d\\mathbf\{x\}\(12\)\+β\|𝒫\|​∫𝒫\(∇⋅𝐮​\(𝐱\)\)2​𝑑𝐱\\displaystyle\+\\frac\{\\beta\}\{\|\\mathcal\{P\}\|\}\\int\_\{\\mathcal\{P\}\}\(\\nabla\\cdot\\mathbf\{u\}\(\\mathbf\{x\}\)\)^\{2\}d\\mathbf\{x\}\(13\)\+γ2​\|Ω\|​∑j∈\{x,y\}∫Ω\(uj​\(𝐱\)−ujT​\(𝐱\)\)2​𝑑𝐱\+δ​\(τv−τu\)2\.\\displaystyle\+\\frac\{\\gamma\}\{2\|\\Omega\|\}\\sum\_\{j\\in\\\{x,y\\\}\}\\int\_\{\\Omega\}\\left\(u\_\{j\}\(\\mathbf\{x\}\)\-u\_\{j\}^\{T\}\(\\mathbf\{x\}\)\\right\)^\{2\}d\\mathbf\{x\}\+\\delta\(\\tau\_\{v\}\-\\tau\_\{u\}\)^\{2\}\.
whereΩ\\Omegadenotes the full computational domain and\|Ω\|\|\\Omega\|its area,ℬ\\mathcal\{B\}denotes the solid obstacle region and\|ℬ\|\|\\mathcal\{B\}\|its area, and𝒫\\mathcal\{P\}denotes the fluid pore space and\|𝒫\|\|\\mathcal\{P\}\|its area\. Furthermore,𝐯=\(vx,vy\)\\mathbf\{v\}=\(v\_\{x\},v\_\{y\}\)is the reference velocity field obtained from the LBM solution,𝐮=\(ux,uy\)\\mathbf\{u\}=\(u\_\{x\},u\_\{y\}\)is the velocity field predicted by the model, andτv\\tau\_\{v\}andτu\\tau\_\{u\}are the reference and predicted tortuosities, respectively\. The coefficientsα\\alpha,β\\beta,γ\\gamma, andδ\\deltaare non\-negative weights controlling the relative contributions of the individual loss components\. All spatial integrals are evaluated numerically on the discrete256×256256\\times 256computational grid\.

The loss function can be decomposed into five contributions\. The first term,ℒvel\\mathcal\{L\}\_\{\\text\{vel\}\}, represents the mean squared error between the predicted and reference velocity fields\. The error is computed separately for each velocity component and then summed\. The second term,ℒobstacle\\mathcal\{L\}\_\{\\text\{obstacle\}\}, imposes an additional penalty when nonzero flow velocity is predicted inside the obstacle regions\. Specifically, within the solid domainℬ\\mathcal\{B\}\(as defined in Equation \([10](https://arxiv.org/html/2605.20250#S2.E10)\)\), the penalty is applied using the L1 norm instead of the L2 norm employed in the velocity loss term\. The use of the L1 norm places greater emphasis on reducing even small nonzero velocity predictions inside obstacles\.

The third term,ℒdiv\\mathcal\{L\}\_\{\\text\{div\}\}, penalizes divergence in the predicted velocity field\. In low\-Mach number flows \(which is our case\), LBM solutions should satisfy the incompressibility condition∇⋅𝐮=0\\nabla\\cdot\\mathbf\{u\}=0everywhere\. This term ensures that the predicted velocity field adheres to this constraint by minimizing divergence errors\.

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/fig_loss_perio.png)Figure 4:Explanation of the loss termℒperio\\mathcal\{L\}\_\{\\text\{perio\}\}defined in equation \([14](https://arxiv.org/html/2605.20250#S2.E14)\)\. a\) The structureSSand translation vectorTT\. b\) Velocity field predicted for structureSS– this velocity field is denoted asC​N​N​\(S\)CNN\(S\)\. c\) Velocity field after imposing translationTT\. d\) StructureSSafter translation by vectorTT: i\.e\.T​\[S\]T\[S\]\. e\) Velocity fieldC​N​N​\(T​\[S\]\)CNN\(T\[S\]\)predicted for translated structureT​\[S\]T\[S\]\. Finallyℒperio\\mathcal\{L\}\_\{\\text\{perio\}\}is computed as integrated squared difference between velocity fields predicted for transformed structureT​\[S\]T\[S\]and transformed velocity field predicted for initial structureSS\.The fourth term,ℒperio\\mathcal\{L\}\_\{\\text\{perio\}\}, addresses the periodicity of the system\. Convolutional neural networks \(CNNs\), by design, are not inherently aware of periodic boundary conditions\. They typically implement border treatment via padding methods—such as zero, replicate, reflect, or circular—which are selected heuristically and do not enforce any physical periodicity\. Several studies demonstrate that the choice of padding can significantly influence prediction accuracy, and that even circular padding fails to generalize across periodic domains without additional architectural constraints or context informationInnamorati et al\. \([2019](https://arxiv.org/html/2605.20250#bib.bib16)\); Gao et al\. \([2021](https://arxiv.org/html/2605.20250#bib.bib7)\); Alguacil et al\. \([2021](https://arxiv.org/html/2605.20250#bib.bib1)\)\. To mitigate this limitation and foster consistency under periodic shifts, we introduce an additional loss termℒperio\\mathcal\{L\}\_\{\\text\{perio\}\}that penalizes discrepancies between prediction for translated structureT​\[S\]T\[S\]and translated prediction for structureSS\. The mechanism by which periodic consistency is promoted through this loss term is illustrated in Figure[4](https://arxiv.org/html/2605.20250#S2.F4)\. The translationTTtakes into account periodic conditions\. For a given structureSSand periodic translationTT, this term is defined as

ℒperio=12​\|Ω\|​∑j∈\{x,y\}∫Ω\(uj​\(𝐱\)−ujT​\(𝐱\)\)2​𝑑𝐱\.\\mathcal\{L\}\_\{\\text\{perio\}\}=\\frac\{1\}\{2\|\\Omega\|\}\\sum\_\{j\\in\\\{x,y\\\}\}\\int\_\{\\Omega\}\\left\(u\_\{j\}\(\\mathbf\{x\}\)\-u\_\{j\}^\{T\}\(\\mathbf\{x\}\)\\right\)^\{2\}d\\mathbf\{x\}\.\(14\)Here,𝐮T\\mathbf\{u\}^\{T\}denotes the prediction obtained after applying a periodic translation to the input structure and translating the resulting predicted velocity field back to the original coordinate system\. This encourages the network to produce predictions that are more consistent across periodic boundaries and better aligned with the periodic assumptions used in the LBM simulations\. During training, several choices of the translation vectorTTwere tested, including random translations\. In the final implementation, we usedT=\[L/2,L/2\]T=\[L/2,L/2\], which provided a practical compromise between training cost and prediction accuracy\.

Finally, the fifth term,ℒtort\\mathcal\{L\}\_\{\\text\{tort\}\}, penalizes discrepancies between the predicted tortuosityτu\\tau\_\{u\}and the reference tortuosityτv\\tau\_\{v\}, which is computed from the ground\-truth LBM solution\. Tortuosity is a key macroscopic transport property that quantifies the effective elongation of flow paths due to the geometry of the porous medium\. Accurate estimation of tortuosity is crucial in many real\-world applications—including filtration, catalysis, and subsurface flow—where it directly impacts performance and efficiency metrics\. By explicitly incorporating tortuosity into the loss function, we ensure that the model not only captures local velocity features but also respects global structural characteristics of the flow field, thereby enhancing the reliability of the model for downstream physical analysis and design tasks\.

All the integrals in Equations \([11](https://arxiv.org/html/2605.20250#S2.E11)\) and \([14](https://arxiv.org/html/2605.20250#S2.E14)\) were approximated by discrete sums over all pixels\. Similarly the divergence inℒdiv\\mathcal\{L\}\_\{\\text\{div\}\}i\.e\.∇⋅𝐮=∂ux∂x\+∂uy∂y\\nabla\\cdot\\mathbf\{u\}=\\frac\{\\partial u\_\{x\}\}\{\\partial x\}\+\\frac\{\\partial u\_\{y\}\}\{\\partial y\}was approximated by first\-order discrete derivatives\.

We emphasize that the model is trained in a supervised manner using velocity fields generated by the LBM solver\. Thus, the network should be interpreted as a surrogate model for the numerical solution under the considered assumptions, rather than as a solver that independently derives the governing equations\. In this work, the term physics\-informed refers to the use of additional physically motivated penalty terms in the loss function, which guide the supervised learning process toward predictions that better satisfy local and global physical constraints\.

### Components of the loss function \- analysis

Table 1:Performance metrics for models trained with different combinations of loss function components\. Each model shares the same ResNet\-101 architecture and dataset\. The table reports the average value of each individual loss component as well as the coefficient of determination \(Rτ2R^\{2\}\_\{\\tau\}\) for tortuosity predictionTable[1](https://arxiv.org/html/2605.20250#S2.T1)presents a quantitative comparison of models trained with different configurations of the custom loss function\. All models share the same ResNet\-101 architecture, which was identified as the most performant among those evaluated, and were trained, validated, and tested on identical datasets\. The values reported represent the average magnitude of each loss component on the test set\. Additionally, the accuracy of tortuosity prediction is evaluated using the coefficient of determinationRτ2R^\{2\}\_\{\\tau\}\. The baseline model was trained using only the velocity loss component, i\.e\.,ℒ=ℒvel\\mathcal\{L\}=\\mathcal\{L\}\_\{\\text\{vel\}\}\. This model achieved the lowest mean squared error \(MSE\) in terms of velocity prediction but failed to suppress spurious flow within obstacle regions, as illustrated in Figure[5](https://arxiv.org/html/2605.20250#S2.F5)\. Furthermore, it demonstrated poor accuracy in predicting tortuosity, withRτ2=0\.127R^\{2\}\_\{\\tau\}=0\.127\.

In the second experiment, we introduced the obstacle loss termℒobstacle\\mathcal\{L\}\_\{\\text\{obstacle\}\}with a weighting factor ofα=5\\alpha=5\. This addition significantly reduced flow predictions within solid regions, with the averageℒobstacle\\mathcal\{L\}\_\{\\text\{obstacle\}\}decreasing by over two orders of magnitude\. It also led to a substantial improvement in tortuosity prediction, increasingRτ2R^\{2\}\_\{\\tau\}to 0\.901 and reducing the MSE inℒtort\\mathcal\{L\}\_\{\\text\{tort\}\}by one order of magnitude\. However, the divergence lossℒdiv\\mathcal\{L\}\_\{\\text\{div\}\}remained similar due to the unconstrained divergence of the predicted field\. To address this, the third model incorporated the divergence loss termℒdiv\\mathcal\{L\}\_\{\\text\{div\}\}withβ=1\\beta=1\. This further reduced the average divergence in the predicted velocity fields and yielded improvements in bothℒvel\\mathcal\{L\}\_\{\\text\{vel\}\}andℒobstacle\\mathcal\{L\}\_\{\\text\{obstacle\}\}, indicating better generalization in void regions and at fluid–solid interfaces\. Despite these improvements, qualitative analysis \(see Figure[5](https://arxiv.org/html/2605.20250#S2.F5)\) revealed that the model occasionally failed to enforce periodic boundary conditions\. To address this issue, the fourth model introduced the periodic consistency lossℒperio\\mathcal\{L\}\_\{\\text\{perio\}\}, weighted byγ=0\.1\\gamma=0\.1\. This modification improved the model’s ability to handle periodic translations, reducingℒperio\\mathcal\{L\}\_\{\\text\{perio\}\}and ensuring greater physical consistency across domain boundaries\. Finally, the fifth and most complete model incorporated the tortuosity consistency lossℒtort\\mathcal\{L\}\_\{\\text\{tort\}\}, with a weighting factor ofδ=0\.01\\delta=0\.01\. This led to the most accurate tortuosity predictions, achieving anRτ2R^\{2\}\_\{\\tau\}value of 0\.983\. Additionally, it exhibited the lowest penalties for both obstacle violation and divergence, while maintaining a velocity prediction error comparable to the baseline\. The final model thus represents the best trade\-off among local accuracy, physical plausibility, and global flow characteristics\.

The above analysis illustrates the contribution of each individual loss function component to the overall model performance\. The sequential inclusion of these terms not only often improves local accuracy in velocity prediction but also enforces key physical constraints such as impermeability of solid regions, incompressibility, and periodic boundary conditions, while enhancing the fidelity of global transport metrics like tortuosity\. Based on these findings, all subsequent models presented in this study were trained using the full loss function with fixed weighting coefficients:α=5\\alpha=5,β=1\\beta=1,γ=0\.1\\gamma=0\.1, andδ=0\.01\\delta=0\.01\.

A more detailed sensitivity analysis of these hyperparameters, including their interactions, is provided in the Supplementary Information \(Section S3\)\.

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/loss_function_components_MM.png)Figure 5:\(left\) Reference velocity field obtained from the LBM simulation\. \(middle\) Prediction from a model trained with a minimal loss function configuration \(α=β=γ=δ=0\\alpha=\\beta=\\gamma=\\delta=0\), illustrating several failure modes\. The red box indicates flow erroneously penetrating solid obstacles\. The green box highlights spurious flow induced by the model’s failure to respect periodic boundary conditions\. The purple box shows a region confined by obstacles on three sides, where the model incorrectly predicts the location of the outflow\. These prediction artifacts are fixed when the full loss function is used, with optimized weights:α=5\\alpha=5,β=1\\beta=1,γ=0\.1\\gamma=0\.1, andδ=0\.01\\delta=0\.01as shown in the right image\.
### Training protocol

During training, the dataset was randomly split into training, validation, and test sets, comprising 70%, 15%, and 15% of the structures, respectively\. The model’s weights were initialized randomly\. We employed a two\-stage training strategy\. In the first stage, a pilot model was trained with divergence, periodicity and tortuosity contributions to total loss set to zero, i\.e\.β=γ=δ=0\\beta=\\gamma=\\delta=0in \([11](https://arxiv.org/html/2605.20250#S2.E11)\)\. In the second stage, the pilot model’s weights were used as a starting point and were then fine\-tuned using the full loss formulation with final values of hyperparameters, i\.e\.β=1,γ=0\.1\\beta=1,\\gamma=0\.1andδ=0\.01\\delta=0\.01\. This staged approach improved training stability and reduced sensitivity to the random dataset split\. To prevent data leakage, both the pilot and final models were trained using the exact same data split\. Optimization was carried out using the Adam algorithmKingma and Ba \([2017](https://arxiv.org/html/2605.20250#bib.bib17)\)\. For the pilot model, a learning rate of 0\.0005 and gradient clippingPascanu et al\. \([2013](https://arxiv.org/html/2605.20250#bib.bib28)\)with a threshold of 1\.0 were used\. For the final model, the learning rate was reduced to10−510^\{\-5\}with gradient clipping at 0\.1\. In both cases, training continued as long as loss measure on validation set decreases \(early stopping technique\)\. Further training details are provided in the Supplementary Information \(Section S2\)\.

Only test set metrics are reported in the paper\. To improve robustness, three different random data splits were evaluated for each model architecture, and the reported results in Table[2](https://arxiv.org/html/2605.20250#Sx1.T2)reflect the average and standard deviation across these splits\.

## Results

Table 2:Performance metrics for various evaluated neural network architectures\. The pixel\-wise root mean squared error \(RMSE\) is reported for velocity field predictions\. For tortuosity and permeability, three metrics are provided: RMSE, mean absolute percentage error \(MAPE\), and coefficient of determination \(R2R^\{2\}\)\. The best result for each metric is highlighted inbold, second\-best withbold\+italic, while the thrid best is shown initalic\. For velocity RMSE standard deviation is reported based on three independent runs with different dataset splits\.Table[2](https://arxiv.org/html/2605.20250#Sx1.T2)presents the predictive performance metrics for all evaluated architectures, covering a range of convolutional backbones with varying levels of architectural complexity and computational efficiency\. The accuracy of the velocity field predictions is assessed using the component\-wise, pixel\-averaged root mean squared error \(RMSE\) defined asRMSE=MSE\\textrm\{RMSE\}=\\sqrt\{\\textrm\{MSE\}\}\. To provide additional context for these values, Figure[2](https://arxiv.org/html/2605.20250#S2.F2)b shows the distributions of both the mean and maximum fluid velocities in the dataset\. The std\(RMSE\) column show the standard deviation of that RMSE computed out of three independent runs with different dataset splits\. The next columns show the performance metrics, i\.e\. the RMSE, Mean Absolute Percentage Error \(MAPE\) and coefficient of determinationR2R^\{2\}for tortuosity and permeability computed from the velocity field predicted by the model and compared with the LBM solution\.

Among all tested models, ResNet\-101 consistently achieves the highest performance across nearly all metrics\. It yields the lowest RMSE for velocity predictions at5\.1⋅10−35\.1\\cdot 10^\{\-3\}, which corresponds to the average magnitude of the prediction error per velocity component\. Given the median fluid velocity of approximately2\.5⋅10−22\.5\\cdot 10^\{\-2\}\(Figure[2](https://arxiv.org/html/2605.20250#S2.F2)b\), this error reflects high predictive accuracy\. Moreover, the low std\(RMSE\) value for ResNet\-101 indicates robustness and stable generalization across dataset variations\.

Architectures designed with a focus on computational efficiency, such as EfficientNet and MobileNetV3, tend to exhibit higher errors across both velocity reconstruction and derived physical quantities\. ConvNeXt, which follows a more modern convolutional design, achieves performance comparable to some of the classical architectures, but does not surpass the best\-performing models\. This trend is consistent across the evaluated metrics, including RMSE, MAPE, andR2R^\{2\}\.

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/lbm_cnn_witherror.png)Figure 6:\(left\) Reference velocity field obtained from LBM simulation\. \(middle\) Prediction from model\. \(right\) Difference in velocity magnitude between reference and predicted velocity fields\. Results for a selected, representative sample for which the RMSE prediction of the velocity field is 5\.8E\-03, minimally greater than the average RMSE of the entire test\-set of the model \(5\.1E\-03, c\.f\. Table[2](https://arxiv.org/html/2605.20250#Sx1.T2)\)\.In Fig\.[6](https://arxiv.org/html/2605.20250#Sx1.F6), we present the prediction results for a selected porous sample in graphical form\. In addition to the velocity field from the LBM method and from the ResNet\-101 prediction, we also show the difference in the absolute value of the velocity field\. A visual analysis of these results confirms that the CNN network can correctly predict fluid flow paths\. Analysis of the error distribution map did not yield any systematic conclusions\. The differences between two methods are not significant and fluctuate around a value of 0, not betraying any trend towards one of the methods\. It is worth noting the excellent continuity of the system at the periodic edges, which is a result of the periodic term in the loss function introduced by us\. Additional error\-bound analysis is provided in the Supplementary Information \(Section S1\)\.

Based on the velocity fields, the tortuosity and permeability of each structure were inferred and compared with the LBM ground truth\. For tortuosity, ResNet\-101 achieves an RMSE of1\.13⋅10−21\.13\\cdot 10^\{\-2\}—small relative to the mean dataset value of 1\.1719 \(Figure[2](https://arxiv.org/html/2605.20250#S2.F2)d\)—and a remarkably low MAPE of 0\.6%\. This indicates that, on average, the predicted tortuosity deviates by less than 1% from the reference\. As shown in Figure[7](https://arxiv.org/html/2605.20250#Sx1.F7)a, the predicted tortuosity values align closely with the exact LBM values, lying along the diagonal with only minor deviations at higher tortuosity values\. The associated error distribution is symmetric and tightly centered around zero\.

For permeability, ResNet\-101 again demonstrates the best performance, achieving the lowest RMSE and MAPE, and the highestR2R^\{2\}across all models\. The MAPE remains below 5%, and the RMSE of8\.69⋅1028\.69\\cdot 10^\{2\}, when compared to the average dataset permeability of2\.3⋅1052\.3\\cdot 10^\{5\}\(Figure[2](https://arxiv.org/html/2605.20250#S2.F2)d\), confirms the accuracy of the predictions\. Figure[7](https://arxiv.org/html/2605.20250#Sx1.F7)b further supports this, showing predicted permeability values tightly clustered around the true LBM values\.

a![Refer to caption](https://arxiv.org/html/2605.20250v1/x2.png)b![Refer to caption](https://arxiv.org/html/2605.20250v1/x3.png)

Figure 7:a\) Tortuosity predicted by the best\-performing model as the function of LBM computed tortuosity\. Inset shows error distribution\. b\) The same for permeability\.Overall, the results in Table[2](https://arxiv.org/html/2605.20250#Sx1.T2)demonstrate that, with the exception of ResNet\-18, most architectures yield comparable performance\. ResNet\-50 consistently performs well and emerges as the second\-best choice, balancing both velocity field prediction and macroscopic property estimation\. The performance achieved by the DenseNet family of convolutional neural networks seems to be systematically lower than that achieved by the simpler ResNet architecture\. Notably, the VGG16 architecture, despite its simpler design and lack of residual or dense connections, achieves competitive results—typically falling between the performance of the ResNet and DenseNet families\. These results suggest that, for the task of mapping pore geometry to velocity fields on a fixed grid, increased architectural complexity or efficiency\-oriented design does not necessarily translate into improved predictive accuracy\.

A comparison with a Fourier Neural OperatorLi et al\. \([2021](https://arxiv.org/html/2605.20250#biba.bib3)\); Kovachki et al\. \([2023](https://arxiv.org/html/2605.20250#biba.bib4)\)model is provided in the Supplementary Information \(Section S4\)\.

### Generalization to different systems

Table 3:Performance metrics of the model on various datasets that are different from the one used to train the model\.While the model is trained on a specific dataset of porous structures, generated using random trigonometric polynomials, the underlying architecture and training methodology can be applied to similar physical systems\. The use of convolutional neural networks allows the model to learn spatially invariant features, enabling it to generalize to geometries and flow patterns that differ from those seen during training\. Moreover, the inclusion of physical constraints in the loss function—such as incompressibility, zero\-flow conditions inside solid regions, and periodicity—encourages the model to learn representations that are physically meaningful rather than dataset\-specific\. It is important to emphasize that all results presented in the following sections were obtained using the best\-performing model, ResNet\-101, trained solely on the original dataset\. No further training, fine\-tuning, or modification was performed for the out\-of\-distribution samples tested here, and the model had no prior exposure to these novel porous structures\. An important question we investigate in the following sections is whether this approach can extend to porous systems with substantially different porosity distributions, obstacle geometries, or boundary conditions, while maintaining predictive accuracy and physical plausibility\. To this end, additional porous structures—distinct from those in the original training dataset—were generated to serve as out\-of\-distribution test cases\. The velocity fields predicted by the model for these new samples were then compared with the corresponding LBM solutions to assess the model’s generalization capabilities\.

#### Obstacles of different shape

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/gen_squares_circles_100_1.png)

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/gen_squares_circles_050_1.png)

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/gen_squares_circles_000_1.png)

Figure 8:LBM computed and neural\-network predicted velocity fields for exemplary systems with circular and square obstacles\. a\) Circle\-only structures \(Pcircle=1\.0P\_\{\\text\{circle\}\}=1\.0\), b\) Mixed structures \(Pcircle=0\.5P\_\{\\text\{circle\}\}=0\.5\), c\) Square\-only structures \(Pcircle=0\.0P\_\{\\text\{circle\}\}=0\.0\)\.A set of porous structures was generated using randomly placed circles and squares\. Obstacle sizes were sampled uniformly between 3 and 8 pixels, and objects were added until the target porosity—randomly chosen between 0\.70 and 0\.95—was achieved\. Each obstacle was a circle with probabilityPcircleP\_\{\\text\{circle\}\}or a square with probability1−Pcircle1\-P\_\{\\text\{circle\}\}\. Three datasets were created withPcircle=0\.00,0\.50,1\.00P\_\{\\text\{circle\}\}=\{0\.00,0\.50,1\.00\}, corresponding to square\-only, mixed, and circle\-only configurations\. Figure[8](https://arxiv.org/html/2605.20250#Sx1.F8)illustrates representative examples of the porous geometries, along with the corresponding LBM solutions and model predictions\. Notably, despite the presence of obstacle shapes that differ from those used during model training, the predictions remain closely aligned with the LBM reference solutions, indicating encouraging generalization to this class of modified geometries\. As summarized in Table[3](https://arxiv.org/html/2605.20250#Sx1.T3), the prediction quality shows a clear dependence on obstacle shape\. Specifically, as the proportion of square obstacles increases, the accuracy of the predicted flow fields deteriorates\. For the circle\-only dataset \(Pcircle=1\.00P\_\{\\text\{circle\}\}=1\.00\), the model achieves performance comparable to that on the original training set, with mean squared errors \(MSEs\) of the same order of magnitude and high coefficients of determination:R2=0\.958R^\{2\}=0\.958for tortuosity andR2=0\.994R^\{2\}=0\.994for permeability\. As the proportion of square obstacles increases, all performance metrics gradually degrade; however, even in the most challenging square\-only case \(Pcircle=0\.00P\_\{\\text\{circle\}\}=0\.00\), the model maintains a mean absolute percentage error \(MAPE\) below 2% for tortuosity and below 6% for permeability, indicating overall robust generalization across varying geometric configurations\.

#### Another boundary conditions

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/gen_border_50.png)

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/gen_border_80.png)

Figure 9:LBM computed and neural\-network predicted velocities fields for exemplary porous systems constrained by two linear obstacles \(pipe\) at porosities a\) 0\.896 and b\) 0\.838\.![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/gen_border_central_1.png)

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/gen_border_central_50.png)

Figure 10:The same as Figure[9](https://arxiv.org/html/2605.20250#Sx1.F9)but with linear obstacle placed in the center of the system at porosities a\) 0\.927, b\) 0\.792\.Periodic boundary conditions were modified by introducing solid vertical obstacles along the horizontal edges of the porous domain, aligned with the direction of fluid flow\. While periodicity is preserved along the horizontal axis, the added boundaries effectively transform the system into an infinite channel or pipe reducing the effect of finite\-size anisotropy and an inclination angle between external body force and resulting streamwise velocityKoza et al\. \([2009](https://arxiv.org/html/2605.20250#bib.bib21)\)\. Importantly, the model was not trained on systems with confined pipe\-like geometries or with any constraint that explicitly prohibited flow across the horizontal edges\. Figure[9](https://arxiv.org/html/2605.20250#Sx1.F9)shows the model’s predictions for three example systems with pipe\-like geometries and periodic boundary conditions\. Although the model was not explicitly trained on such configurations, it correctly reconstructs the overall flow profile\. Importantly, the velocity component perpendicular to the pipe axis remains negligible, demonstrating that the model implicitly respects the physical constraint of unidirectional flow in these settings\.

An equivalent physical setup can be achieved by placing a single solid barrier horizontally across the center of the sample\. Although this configuration is physically analogous to side\-wall confinement, it may pose a greater challenge to the model, as it relies more strongly on accurately capturing internal boundary effects rather than edge constraints alone – some exemplary predictions made by the model are shown in Figure[10](https://arxiv.org/html/2605.20250#Sx1.F10)\.

As shown by the performance metrics in Table[3](https://arxiv.org/html/2605.20250#Sx1.T3), the model maintains high predictive accuracy for the in\-pipe systems, with only a slight reduction in performance compared to the original dataset\. The coefficient of determination \(R2R^\{2\}\) exceeds 0\.96 for tortuosity and 0\.99 for permeability, while the mean absolute percentage error \(MAPE\) remains below 0\.8% and 4\.5%, respectively\. These results indicate that the model generalizes well to previously unseen pipe\-like geometries and continues to provide accurate predictions for both macroscopic transport properties in these test cases\.

### Generalization to different porosities

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/gen_porosity07_1.png)

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/gen_porosity06_1.png)

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/gen_porosity05_1.png)

Figure 11:LBM computed and neural\-network predicted velocities fields for exemplary porous systems of lower porosity\. Systems have porosities: a\) 0\.752, b\) 0\.649, c\) 0\.512\.Three additional datasets were constructed to evaluate the model’s performance on structures with decreasing porosity\. The first dataset includes samples with porosity in the range0\.7<φ<0\.80\.7<\\varphi<0\.8, the second with0\.6<φ<0\.70\.6<\\varphi<0\.7, and the third with0\.5<φ<0\.60\.5<\\varphi<0\.6\. Each dataset consists of 500 structures generated using the method defined in Equation \([9](https://arxiv.org/html/2605.20250#S2.E9)\), with the threshold parameterε\\varepsilonin Equation \([10](https://arxiv.org/html/2605.20250#S2.E10)\) adjusted to achieve the desired porosity\. Figure[11](https://arxiv.org/html/2605.20250#Sx1.F11)shows representative LBM solutions and corresponding predictions made by the model, each taken from a randomly selected sample within the respective datasets\.

Performance metrics for each case are summarized in Table[3](https://arxiv.org/html/2605.20250#Sx1.T3)\. The results reveal a clear dependency of model accuracy on porosity\. As porosity decreases, prediction performance deteriorates—as reflected by decliningR2R^\{2\}scores and increasing MAPE for both tortuosity and permeability\. Notably, the first dataset \(0\.7<φ<0\.80\.7\\\!<\\\!\\varphi\\\!<\\\!0\.8\) falls within the porosity range used during model training, and the model performs well on these samples\. For the second dataset \(0\.6<φ<0\.70\.6\\\!<\\\!\\varphi\\\!<\\\!0\.7\), which lies outside the training distribution, the model still achieves reasonably high accuracy, withR2≈0\.8R^\{2\}\\approx 0\.8for tortuosity and over 0\.9 for permeability, while MAPE values remain within acceptable bounds\.

However, for the most dense structures \(0\.5<φ<0\.60\.5\\\!<\\\!\\varphi\\\!<\\\!0\.6\), a substantial drop in performance is observed across all metrics\. As the system approaches the percolation threshold, where critical phenomena begin to dominate transport behavior, the model’s predictive capability breaks down\. This highlights a fundamental limitation in the current approach—its reduced transferability in regimes exhibiting extreme porosity and critical connectivity transitions\.

The purpose of this experiment was to assess out\-of\-distribution generalization without modifying the training set, rather than to claim reliable performance across all porosity regimes\. Low\-porosity systems are particularly challenging because transport becomes increasingly controlled by narrow channels, connectivity effects, and proximity to the percolation threshold\. The observed deterioration therefore identifies a limitation of the present model\. Improving accuracy in this regime would likely require enriching the training data with additional low\-porosity samples, for example through adaptive sampling or curriculum learning strategies\.

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/lio2_806.png)

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/lio2_745.png)

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/lio2_796.png)

Figure 12:LBM computed and neural\-network predicted velocities fields for exemplary real Li\-O2electrodes\. Systems have porosities: a\) 0\.728, b\) 0\.771, c\) 0\.802\.
### Real porous structures

To evaluate the model on highly irregular real geometries, we tested it on the Li\-O2electrodes dataset \(DRP\-1129v2 from theDigital Porous Media Portal \([2024](https://arxiv.org/html/2605.20250#bib.bib4)\); Prodanovic et al\. \([2025](https://arxiv.org/html/2605.20250#bib.bib29)\); Turhan et al\. \([2024](https://arxiv.org/html/2605.20250#bib.bib35)\)\)\. This dataset consists of three\-dimensional lithium\-ion battery electrode microstructures reconstructed from Nano\-CT scans, with porosity in the range of approximately 0\.73–0\.91\. These structures differ substantially from the synthetic training data, exhibiting complex, heterogeneous pore networks formed through physical manufacturing processes\. In our study, two\-dimensional slices were extracted from the volumetric data and used as inputs, with corresponding velocity fields computed using LBM simulations\.

The results, summarized in Table[3](https://arxiv.org/html/2605.20250#Sx1.T3), show that the model maintains reasonable predictive accuracy despite the significant domain shift, with a velocity RMSE of6\.41⋅10−36\.41\\cdot 10^\{\-3\}and high agreement for permeability \(R2=0\.986R^\{2\}=0\.986\)\. At the same time, the accuracy for tortuosity decreases \(R2=0\.798R^\{2\}=0\.798\), reflecting the increased geometric complexity\. These findings indicate that the model generalizes to realistic, highly irregular porous structures, while also highlighting the expected degradation in performance outside the training distribution\. Figure[12](https://arxiv.org/html/2605.20250#Sx1.F12)shows representative LBM solutions and corresponding predictions made by the model\.

### 2\.1Use of CNN predictions for LBM initialization

One of the problems with calculating flows in porous media using the LBM method at extremely high and low porosities is the hindered convergence of the method due to the diffusive nature of momentum transportMatyka et al\. \([2008](https://arxiv.org/html/2605.20250#bib.bib26)\)\. The classic, so\-called ’cold’ start of the LBM method involves running the calculation from a zero\-velocity field\. It turns out that even the imperfect velocity field predictions from neural networks can be an effective remedy to this problem and serve as a starting condition for the LBM method, which adds to the system at the end, smoothing out errors\.

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/speedup.png)Figure 13:The graph presents the number of steps of the LBM when starting from the output of the CNN \(so\-called warm start, x axis\) to the number of steps whether the LBM starts from zero \(so\-called cold start, y axis\)\. All points below the line mean the LBM\+CNN procedure gain time\. Points are coloured by porosity\.We decided to test this method on samples outside the scope of network learning and perform a test in which a neural network trained on trigonometric data makes predictions of flows in arrangements of overlapping figures, and then feeds such prepared data back to the LBM solver\.

As a result of starting from an approximate velocity field generated by the CNNs, the LBM method in90%90\\%systems converged faster to the correct solution \(see Figure[13](https://arxiv.org/html/2605.20250#Sx1.F13)\), compared to the solver that started from a zero velocity field\. The median reduction in the number of iterations was approximately50%50\\%, meaning that in half of the tested systems the reduction was greater than50%50\\%, and in half it was smaller\. The interquartile range, given byQ1≈22%Q\_\{1\}\\approx 22\\%andQ3≈68%Q\_\{3\}\\approx 68\\%, indicates that50%50\\%of the systems exhibit reductions within this interval\. A bootstrap procedure \(based on repeated resampling of the dataset and recomputation of the median\) yields a95%95\\%confidence interval of approximately\[46%,54%\]\[46\\%,54\\%\], which quantifies the statistical uncertainty of the estimated median reduction and shows that this value is stable across the dataset\. A paired Wilcoxon signed\-rank test confirms that this improvement is highly statistically significant \(p≈4×10−70p\\approx 4\\times 10^\{\-70\}\)\.

With the additional observation that the trained neural network produces results orders of magnitude faster, these findings provide hope for the practical application of our results\.

### 2\.2Computational time comparison

Neural network enables fast inference\. The average prediction time is just 5 ms per structure on GPU \(NVIDIA GeForce RTX 4090\), while on a CPU \(AMD Ryzen 9 5950X\) it averages 190 ms\. Importantly, the inference time remains consistent regardless of porosity or other structural properties\. In comparison, performing LBM calculations the same set of structures requires an average of 560 ms per structure \(utilizing the same CPU and the same number of cores\)\. It should be noted that this comparison concerns inference only\. Training is computationally expensive but performed once and amortized over repeated predictions\.

## 3Conclusions

We have presented a complete pipeline for predicting pore\-scale fluid flow through porous media using physics\-informed convolutional neural networks\. The main contribution of this work is the construction and evaluation of a custom loss function tailored to this problem, combining velocity\-field reconstruction with penalties for flow inside solid regions, incompressibility, periodic consistency, and tortuosity matching\. Our analysis of weights for loss function components showed complex dependency of each part on the overall error of fluid flow predictions with the winning combination ofα=5,β=1,γ=0\.1\\alpha=5,\\beta=1,\\gamma=0\.1andδ=0\.01\\delta=0\.01\. We suggest there is still space to optimize both the choice of the neural network architecture \(ResNet\-101 in our case\), its size and combination of loss weights\. Future work may further improve performance through systematic architecture search, automated loss\-weight optimization, or adaptive sampling of challenging regimes\.

The key finding of this work, however, is generalization\. While the model is trained on a specific dataset of porous structures, the underlying architecture and training methodology are designed to be broadly applicable to a wide class of similar physical systems\. The use of convolutional neural networks allows the model to learn spatially invariant features, enabling it to generalize to geometries and flow patterns that differ from those seen during training\. Moreover, the inclusion of physical constraints in the loss function\-such as incompressibility, zero\-flow conditions inside solid regions, and periodicity\-encourages the model to learn representations that are physically meaningful rather than overly dataset\-specific\. As a result, the trained model can be adapted to other porous systems, including those with different porosity distributions, obstacle shapes, or boundary conditions, provided that they share the same fundamental flow characteristics\. For significantly different domains, transfer learning or fine\-tuning on a small number of new samples may be sufficient to adapt the model while retaining much of the learned structure\. Future work could explore this direction in more detail by systematically evaluating the model’s performance on out\-of\-distribution datasets or synthetic benchmarks with known analytical properties\.

Despite the encouraging results, several limitations remain\. The present implementation is restricted to two\-dimensional, fixed\-resolution geometries, and its accuracy decreases for strongly out\-of\-distribution cases, such as low\-porosity structures near connectivity transitions\. However, the approach is, by design, naturally extendable to three\-dimensional systems: the same encoder–decoder concept, backbone\-based architecture, and physics\-informed loss strategy can be retained, with two\-dimensional convolutional operations replaced by their three\-dimensional counterparts\. Moreover, the physical constraints are imposed through penalty terms and are therefore only approximately enforced\. Future work should explore full three\-dimensional implementations, transfer learning for experimentally obtained microstructures, systematic optimization of loss weights, and possible integration with other physics\-informed or operator\-learning frameworks, such as PINNs, Fourier Neural Operators, or DeepONets\.

## 4Code and data availability

Both the dataset and the code are publicly available in a dedicated[GitHub repository](https://github.com/dioscuri-tda/PhysInfPorousFluidFlow)\. Moreover, we provide the best\-performing pretrained model, which can be used directly without requiring retraining\.

## Author Contributions

R\.T\., P\.D\. and M\.M\. started and designed the project\. R\.T\. developed codebase for neural networks, train the models and analyzed the results M\.M\. prepared fluid flow simulations codes\. R\.T\. and M\.M\. designed the loss function to incorporate flow conditions into training\. P\.D\. suggested time analysis and LBM initialization procedure\. R\.T\. and M\.M\. wrote the first draft of manuscript\. All authors participated in writing the text and discussing the results\.

## 5Funding

R\. Topolnicki and P\. Dłotko acknowledges the support of Dioscuri program initiated by the Max Planck Society, jointly managed with the National Science Centre \(Poland\), and mutually funded by the Polish Ministry of Science and Higher Education and the German Federal Ministry of Education and Research\. Funded by National Science Centre, Poland under the OPUS call in the Weave programme 2021/43/I/ST3/00228\. This research was funded in whole or in part by National Science Centre \(2021/43/I/ST3/00228\)\. For the purpose of Open Access, the author has applied a CC\-BY public copyright licence to any Author Accepted Manuscript \(AAM\) version arising from this submission\.

M\. Matyka acknowledges the financial support from the Slovenian Research And Innovation Agency \(ARIS\) research core funding No\. P2\-0095\.

Authors acknowledge the financial support from M\-ERA\.NET’s PORMETALOMICS project supported by MCIN/AEI/10\.13039/501100011033 and the European Union’s NextGenerationEU/PRTR funds\. Financial support from the PORMETALOMICS project, funded by the National Science Centre, Poland \(project no\. 2021/03/Y/ST5/00232\) within the M\-ERA\.NET 3 programme, is also gratefully acknowledged\. This project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 958174\.

## Acknowledgements

The computations in this work were partially supported by the Center for Artificial Intelligence at Adam Mickiewicz University\. We would like to thank Prof\. K\. Jassem and Dr\. B\. Naskrecki for granting us access to the Center’s infrastructure\. The calculations were made with the support of the Interdisciplinary Center for Mathematical and Computational Modeling of the University of Warsaw \(ICM UW\) under the computational grant no g105\-2705\.

## References

- Alguacil et al\. \(2021\)Alguacil, A\., Pinto, W\.G\., Bauerheim, M\., Jacob, M\.C\., Moreau, S\., 2021\.Effects of boundary conditions in fully convolutional networks for learning spatio\-temporal dynamics\.URL:[https://arxiv\.org/abs/2106\.11160](https://arxiv.org/abs/2106.11160),[arXiv:2106\.11160](http://arxiv.org/abs/2106.11160)\.
- Anderson \(2002\)Anderson, J\.D\., 2002\.Computational fluid dynamics: the basics with applications\.McGraw\-Hill New York\.
- Cossins et al\. \(2023\)Cossins, T\., Mishra, A\., Haese, R\.R\., 2023\.The feasibility of enhanced pore space utilization in co2 storage reservoirs using an artificially emplaced si\-gel flow barrier\.Scientific Reports 13, 9334\.
- Digital Porous Media Portal \(2024\)Digital Porous Media Portal, 2024\.DRP\-1129v2 dataset\.URL:[https://doi\.org/10\.17612/dpm/DRP\-1129v2](https://doi.org/10.17612/dpm/DRP-1129v2),[10\.17612/dpm/DRP\-1129v2](http://dx.doi.org/10.17612/dpm/DRP-1129v2)\.
- Duda et al\. \(2011\)Duda, A\., Koza, Z\., Matyka, M\., 2011\.Hydraulic tortuosity in arbitrary porous media flow\.Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 84, 036319\.
- Fong \(2015\)Fong, C\.W\., 2015\.Permeability of the blood–brain barrier: molecular mechanism of transport of drugs and physiologically important compounds\.The Journal of membrane biology 248, 651–669\.
- Gao et al\. \(2021\)Gao, H\., Sun, L\., Wang, J\.X\., 2021\.Phygeonet: Physics\-informed geometry\-adaptive convolutional neural networks for solving parameterized steady\-state pdes on irregular domain\.Journal of Computational Physics 428, 110079\.doi:[https://doi\.org/10\.1016/j\.jcp\.2020\.110079](http://dx.doi.org/https://doi.org/10.1016/j.jcp.2020.110079)\.
- Graczyk and Matyka \(2020\)Graczyk, K\.M\., Matyka, M\., 2020\.Predicting porosity, permeability, and tortuosity of porous media from images by deep learning\.Scientific reports 10, 21488\.
- Graczyk et al\. \(2023\)Graczyk, K\.M\., Strzelczyk, D\., Matyka, M\., 2023\.Deep learning for diffusion in porous media\.Scientific Reports 13, 9769\.
- Guo et al\. \(2024\)Guo, C\., Wang, Y\., Han, Y\., Ji, M\., Wu, Y\., 2024\.Unsteady flow\-field forecasting leveraging a hybrid deep\-learning architecture\.Physics of Fluids 36\.
- Guo et al\. \(2016\)Guo, X\., Li, W\., Iorio, F\., 2016\.Convolutional neural networks for steady flow approximation, in: Proceedings of the 22nd ACM SIGKDD international conference on knowledge discovery and data mining, pp\. 481–490\.
- Guo and Shu \(2013\)Guo, Z\., Shu, C\., 2013\.Lattice Boltzmann method and its application in engineering\. volume 3\.World Scientific\.
- He et al\. \(2015\)He, K\., Zhang, X\., Ren, S\., Sun, J\., 2015\.Deep residual learning for image recognition\.URL:[https://arxiv\.org/abs/1512\.03385](https://arxiv.org/abs/1512.03385),[arXiv:1512\.03385](http://arxiv.org/abs/1512.03385)\.
- Howard et al\. \(2019\)Howard, A\., Sandler, M\., Chu, G\., Chen, L\.\-C\., Chen, B\., Tan, M\., Wang, W\., Zhu, Y\., Pang, R\., Vasudevan, V\., Le, Q\.V\., Adam, H\., 2019\.Searching for MobileNetV3\.URL:[https://arxiv\.org/abs/1905\.02244](https://arxiv.org/abs/1905.02244),[arXiv:1905\.02244](http://arxiv.org/abs/1905.02244)\.
- Huang et al\. \(2018\)Huang, G\., Liu, Z\., van der Maaten, L\., Weinberger, K\.Q\., 2018\.Densely connected convolutional networks\.URL:[https://arxiv\.org/abs/1608\.06993](https://arxiv.org/abs/1608.06993),[arXiv:1608\.06993](http://arxiv.org/abs/1608.06993)\.
- Innamorati et al\. \(2019\)Innamorati, C\., Ritschel, T\., Weyrich, T\., Mitra, N\.J\., 2019\.Learning on the edge: Investigating boundary filters in cnns\.International Journal of Computer Vision 128, 773–782\.doi:[10\.1007/s11263\-019\-01223\-y](http://dx.doi.org/10.1007/s11263-019-01223-y)\.
- Kingma and Ba \(2017\)Kingma, D\.P\., Ba, J\., 2017\.Adam: A method for stochastic optimization\.URL:[https://arxiv\.org/abs/1412\.6980](https://arxiv.org/abs/1412.6980),[arXiv:1412\.6980](http://arxiv.org/abs/1412.6980)\.
- Koponen et al\. \(1996\)Koponen, A\., Kataja, M\., Timonen, J\., 1996\.Tortuous flow in porous media\.Physical Review E 54, 406\.
- Koponen et al\. \(1997\)Koponen, A\., Kataja, M\., Timonen, J\., 1997\.Permeability and effective porosity of porous media\.Physical Review E 56, 3319\.
- Kovachki et al\. \(2023\)Kovachki, N\., Li, Z\., Azizzadenesheli, K\., Liu, B\., Bhattacharya, K\., Stuart, A\., Anandkumar, A\., 2023\.Neural Operator: Learning Maps Between Function Spaces\.Journal of Machine Learning Research 24, 1–97\.
- Koza et al\. \(2009\)Koza, Z\., Matyka, M\., Khalili, A\., 2009\.Finite\-size anisotropy in statistically uniform porous media\.Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 79, 066306\.
- Lehmann et al\. \(2022\)Lehmann, M\., Krause, M\.J\., Amati, G\., Sega, M\., Harting, J\., Gekle, S\., 2022\.Accuracy and performance of the lattice boltzmann method with 64\-bit, 32\-bit, and customized 16\-bit number formats\.Physical Review E 106, 015308\.
- Li et al\. \(2021\)Li, Z\., Kovachki, N\., Azizzadenesheli, K\., Liu, B\., Bhattacharya, K\., Stuart, A\., Anandkumar, A\., 2021\.Fourier Neural Operator for Parametric Partial Differential Equations\.URL:[https://arxiv\.org/abs/2010\.08895](https://arxiv.org/abs/2010.08895),[arXiv:2010\.08895](http://arxiv.org/abs/2010.08895)\.
- Lin et al\. \(2025\)Lin, Y\., Wu, Z\., You, S\., Yang, C\., Wang, Q\., Yin, W\., Qiu, T\., 2025\.Deep learning\-based prediction of velocity and temperature distributions in metal foam with hierarchical pore structure\.Green Chemical Engineering 6, 209–222\.
- Liu et al\. \(2022\)Liu, Z\., Mao, H\., Wu, C\.\-Y\., Feichtenhofer, C\., Darrell, T\., Xie, S\., 2022\.A ConvNet for the 2020s\.URL:[https://arxiv\.org/abs/2201\.03545](https://arxiv.org/abs/2201.03545),[arXiv:2201\.03545](http://arxiv.org/abs/2201.03545)\.
- Matyka et al\. \(2008\)Matyka, M\., Khalili, A\., Koza, Z\., 2008\.Tortuosity\-porosity relation in porous media flow\.Physical Review E—Statistical, Nonlinear, and Soft Matter Physics 78, 026306\.
- Matyka and Koza \(2012\)Matyka, M\., Koza, Z\., 2012\.How to calculate tortuosity easily?, in: AIP Conference Proceedings 4, American Institute of Physics\. pp\. 17–22\.
- Pascanu et al\. \(2013\)Pascanu, R\., Mikolov, T\., Bengio, Y\., 2013\.On the difficulty of training recurrent neural networks\.URL:[https://arxiv\.org/abs/1211\.5063](https://arxiv.org/abs/1211.5063),[arXiv:1211\.5063](http://arxiv.org/abs/1211.5063)\.
- Prodanovic et al\. \(2025\)Prodanovic, M\., Esteva, M\., Ketcham, R\., Chang, B\., Turhan, C\., Gentle, J\., Khan, S\., Belcher, V\., 2025\.Digital Porous Media Portal \(DPMP\) for publication, analysis, and simulation of porous media images\.Digital Porous Media Portal\.doi:[10\.17612/FGMN\-D889](http://dx.doi.org/10.17612/FGMN-D889)\.
- Ronneberger et al\. \(2015\)Ronneberger, O\., Fischer, P\., Brox, T\., 2015\.U\-net: Convolutional networks for biomedical image segmentation\.URL:[https://arxiv\.org/abs/1505\.04597](https://arxiv.org/abs/1505.04597),[arXiv:1505\.04597](http://arxiv.org/abs/1505.04597)\.
- Simonyan and Zisserman \(2015\)Simonyan, K\., Zisserman, A\., 2015\.Very deep convolutional networks for large\-scale image recognition\.URL:[https://arxiv\.org/abs/1409\.1556](https://arxiv.org/abs/1409.1556),[arXiv:1409\.1556](http://arxiv.org/abs/1409.1556)\.
- Succi \(2001\)Succi, S\., 2001\.The lattice Boltzmann equation: for fluid dynamics and beyond\.Oxford university press\.
- Tan and Le \(2019\)Tan, M\., Le, Q\.V\., 2019\.EfficientNet: Rethinking Model Scaling for Convolutional Neural Networks, in: Proceedings of the 36th International Conference on Machine Learning\.6105–6114\.
- Tomczak and Szafran \(2019\)Tomczak, T\., Szafran, R\.G\., 2019\.A new gpu implementation for lattice\-boltzmann simulations on sparse geometries\.Computer Physics Communications 235, 258–278\.
- Turhan et al\. \(2024\)Turhan, Ç\., Chang, B\., Mohamed, A\., Esteva, M\., Ketcham, R\., McClure, J\., Prodanovic, M\., 2024\.Digital Porous Media Portal for image curation, characterization, visualization, and transport simulation in porous media, in: International Symposium of the Society of Core Analysts\.SCA2024\-1056\.
- Zhao et al\. \(2023\)Zhao, X\., Gong, Z\., Zhang, Y\., Yao, W\., Chen, X\., 2023\.Physics\-informed convolutional neural networks for temperature field prediction of heat source layout without labeled data\.Engineering Applications of Artificial Intelligence 117, 105516\.
- Zhong et al\. \(2021\)Zhong, L\., Han, X\., Yuan, X\., Liu, Y\., Zou, J\., Wang, Q\., 2021\.Permeability variation and its impact on oil recovery from unconsolidated sand heavy\-oil reservoirs during steamflooding process\.SPE Reservoir Evaluation & Engineering 24, 159–173\.

## Supplementary Information

Supplementary Information Physics\-informed convolutional neural networks for fluid flow through porous media

###### Contents

1. [1Introduction](https://arxiv.org/html/2605.20250#S1)
2. [2Materials and Methods](https://arxiv.org/html/2605.20250#S2)
3. [2\.1Use of CNN predictions for LBM initialization](https://arxiv.org/html/2605.20250#Sx1.SS1)
4. [2\.2Computational time comparison](https://arxiv.org/html/2605.20250#Sx1.SS2)
5. [3Conclusions](https://arxiv.org/html/2605.20250#S3)
6. [4Code and data availability](https://arxiv.org/html/2605.20250#S4)
7. [5Funding](https://arxiv.org/html/2605.20250#S5)
8. [References](https://arxiv.org/html/2605.20250#bib)
9. [Supplementary Information](https://arxiv.org/html/2605.20250#Ax1)
10. [S1Error Bounds](https://arxiv.org/html/2605.20250#A1)
11. [S2Training procedure details](https://arxiv.org/html/2605.20250#A2)
12. [S3Loss function hyperparameters](https://arxiv.org/html/2605.20250#A3)
13. [S4Fourier Neural Operators](https://arxiv.org/html/2605.20250#A4)
14. [References](https://arxiv.org/html/2605.20250#biba)

## Appendix S1Error Bounds

To assess the robustness of the surrogate in the regimes highlighted by the referee, we complement global error metrics with a conditional uncertainty analysis based on the velocity magnitude\|v\|\|v\|\. In Fig\.[S1](https://arxiv.org/html/2605.20250#A1.F1), we plot a two\-dimensional histogram of\|v\|pred\|v\|\_\{\\mathrm\{pred\}\}versus\|v\|LBM\|v\|\_\{\\mathrm\{LBM\}\}, which shows both the concentration of samples around the ideal relationy=xy=xand the spread of the prediction error as a function of the predicted velocity magnitude\. The purpose of estimating confidence bounds is to move beyond a single averaged score and quantify how predictive uncertainty changes across the state space, in particular in difficult low\-porosity regimes where transport pathways become constricted and the flow becomes more heterogeneous\. To this end, we define residuals

r=\|v\|LBM−\|v\|pred,r=\|v\|\_\{\\mathrm\{LBM\}\}\-\|v\|\_\{\\mathrm\{pred\}\},\(S1\)and estimate conditional lower and upper quantiles ofrrwithin bins of the predicted velocity magnitude\|v\|pred\|v\|\_\{\\mathrm\{pred\}\}\. Denoting the conditional quantile function byQα​\(⋅\)Q\_\{\\alpha\}\(\\cdot\), the empirical estimates correspond to

qα\(x\)=Qα\(r\|\|v\|pred=x\),q\_\{\\alpha\}\(x\)=Q\_\{\\alpha\}\\big\(r\\,\\big\|\\,\|v\|\_\{\\mathrm\{pred\}\}=x\\big\),\(S2\)whereα∈\{0\.05,0\.95\}\\alpha\\in\\\{0\.05,0\.95\\\}in the present work\. These quantiles define lower and upper bounds for the ground\-truth velocity magnitude conditioned on the prediction:

v¯​\(x\)=x\+q0\.05​\(x\),v¯​\(x\)=x\+q0\.95​\(x\)\.\\underline\{v\}\(x\)=x\+q\_\{0\.05\}\(x\),\\qquad\\overline\{v\}\(x\)=x\+q\_\{0\.95\}\(x\)\.\(S3\)
The discrete estimates\{\(xi,qα​\(xi\)\)\}\\\{\(x\_\{i\},q\_\{\\alpha\}\(x\_\{i\}\)\)\\\}obtained from binning are then converted into continuous functions using a shape\-preserving piecewise cubic Hermite interpolant \(PCHIP\),Fritsch and Carlson \([1980](https://arxiv.org/html/2605.20250#biba.bib1)\); Fritsch and Butland \([1984](https://arxiv.org/html/2605.20250#biba.bib2)\)\. In practice, this corresponds to constructing a piecewise cubic functionsα​\(x\)s\_\{\\alpha\}\(x\)such that

sα​\(xi\)=qα​\(xi\),i=1,…,N,s\_\{\\alpha\}\(x\_\{i\}\)=q\_\{\\alpha\}\(x\_\{i\}\),\\quad i=1,\\dots,N,\(S4\)with continuity of bothsα​\(x\)s\_\{\\alpha\}\(x\)and its first derivative, while preserving the local monotonicity of the data\. The final confidence bounds are therefore given by

v¯​\(x\)=x\+s0\.05​\(x\),v¯​\(x\)=x\+s0\.95​\(x\)\.\\underline\{v\}\(x\)=x\+s\_\{0\.05\}\(x\),\\qquad\\overline\{v\}\(x\)=x\+s\_\{0\.95\}\(x\)\.\(S5\)
These discrete quantile estimates are then converted into continuous lower and upper bound functions using a shape\-preserving piecewise cubic Hermite interpolant, PCHIP,Fritsch and Carlson \([1980](https://arxiv.org/html/2605.20250#biba.bib1)\); Fritsch and Butland \([1984](https://arxiv.org/html/2605.20250#biba.bib2)\)\. In practice, this interpolation constructs a cubic polynomial within each bin interval while enforcing continuity of both the function and its first derivative across bin boundaries\. Unlike standard cubic splines, which may introduce non\-physical oscillations or overshoots in regions with sparse or rapidly changing data, PCHIP explicitly preserves the local monotonicity and shape of the underlying quantile estimates\. As a result, the interpolated bounds remain consistent with the empirical trends observed in the binned data, providing smooth yet faithful approximations of the conditional quantile curves without introducing artificial extrema or distortions\.

![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/flow_ci_histogram.png)Figure S1:Two\-dimensional histogram of ground\-truth \(LBM\-computed\) velocity magnitude\|v\|LBM\|v\|\_\{\\mathrm\{LBM\}\}versus predicted velocity magnitude\|v\|pred\|v\|\_\{\\mathrm\{pred\}\}\. The color scale represents the density of samples\. The dashed line indicates the ideal relation\|v\|LBM=\|v\|pred\|v\|\_\{\\mathrm\{LBM\}\}=\|v\|\_\{\\mathrm\{pred\}\}, while the shaded region denotes the estimated90%90\\%conditional confidence band for\|v\|LBM\|v\|\_\{\\mathrm\{LBM\}\}given\|v\|pred\|v\|\_\{\\mathrm\{pred\}\}, obtained from empirical quantiles of the residualsr=\|v\|LBM−\|v\|predr=\|v\|\_\{\\mathrm\{LBM\}\}\-\|v\|\_\{\\mathrm\{pred\}\}\. The narrow band in the high\-density region indicates strong agreement between predictions and ground truth for moderate velocities, whereas the widening of the band at higher velocities reflects increased uncertainty in more challenging flow regimes\.The resulting bounds are shown globally as a band around the diagonal in Fig\.[S1](https://arxiv.org/html/2605.20250#A1.F1), and spatially in Fig\.[S2](https://arxiv.org/html/2605.20250#A1.F2), where the lower and upper confidence limits are evaluated pointwise as functions of the predicted velocity magnitude\|v\|pred\|v\|\_\{\\mathrm\{pred\}\}\. The multi\-panel visualization provides additional insight into the spatial structure of the uncertainty: regions of high velocity, typically corresponding to preferential flow channels, remain tightly constrained, whereas low\-velocity regions near solid boundaries and in poorly connected pore spaces exhibit a visibly wider confidence envelope\. These regions are associated with reduced porosity and proximity to the percolation threshold, where the flow becomes more sensitive to small geometric variations\. The predicted fields remain consistent with the topology and spatial organization of the LBM solution, and the ground\-truth solution is, in most regions, contained within the estimated bounds\. This indicates that the model retains meaningful predictive capability while appropriately reflecting increased uncertainty in more complex flow configurations\.

Overall, this conditional quantile\-based analysis provides a more detailed characterization of model performance, highlighting how predictive reliability varies across the domain and offering spatially resolved bounds on the expected prediction error, rather than relying solely on global summary metrics\.

a\)![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/flow_ci_205.png) b\)![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/flow_ci_100.png) c\)![Refer to caption](https://arxiv.org/html/2605.20250v1/figures/flow_ci_101.png)

Figure S2:Spatial visualization of conditional confidence bounds for the velocity magnitude\|v\|\|v\|\. The left panel shows the ground\-truth LBM solution, while the central and right panels display the estimated lower and upper confidence bounds, respectively, obtained from the conditional quantile model\. Solid regions correspond to the pore structure\. The bounds are evaluated pointwise as functions of the predicted velocity magnitude\|v\|pred\|v\|\_\{\\mathrm\{pred\}\}, resulting in spatially varying uncertainty estimates\. Systems have porosities: a\) 0\.713, b\) 0\.817, c\) 0\.931\.The same uncertainty quantification framework can be applied independently to the individual components of the velocity vector\. In this case, the analysis is performed for each component \(e\.g\.,vxv\_\{x\}andvyv\_\{y\}\) by constructing conditional confidence intervals for the ground\-truth component given its predicted value\. This involves defining residualsr=vLBM−vpredr=v\_\{\\mathrm\{LBM\}\}\-v\_\{\\mathrm\{pred\}\}for each component and estimating conditional quantiles ofrrwithin bins of the predicted component\. This extension allows for component\-wise uncertainty characterization, capturing potential anisotropies in the prediction error while preserving the local structure of the flow field\. In practice, the results are qualitatively consistent with the magnitude\-based analysis, while providing additional directional insight into the model performance\.

## Appendix S2Training procedure details

Training was performed on a workstation equipped with an NVIDIA GeForce RTX 4090 GPU with 24 GB of VRAM and an AMD Ryzen 9 5950X CPU\. The dataset was randomly split into training, validation, and test subsets containing 70%, 15%, and 15% of the samples, respectively\. A batch size of 32 was used in all experiments\. During training, the loss was monitored on the validation set, and early stopping was applied with a patience of 20 epochs\. The total number of epochs varied between architectures, but training typically required approximately 250 epochs\. Optimization was performed using the Adam algorithm with an initial learning rate of5×10−45\\times 10^\{\-4\}\. The learning rate was updated using a step scheduler with step size 10 and decay factorgamma=0\.8\.

We found that the most stable convergence was obtained using a two\-stage training procedure\. In the first stage, the model was trained with the obstacle loss enabled and the remaining physics\-based terms disabled, i\.e\. withα=5\\alpha=5andβ=γ=δ=0\\beta=\\gamma=\\delta=0\. In the second stage, training was resumed from the weights obtained in the first stage using the full loss function with the final hyperparameter valuesα=5,β=1,γ=0\.1\\alpha=5,\\beta=1,\\gamma=0\.1andδ=0\.01\\delta=0\.01\.

## Appendix S3Loss function hyperparameters

The loss function introduced in Eq\. \(11\) of the manuscript contains four hyperparameters:α\\alpha,β\\beta,γ\\gamma, andδ\\delta, which control the relative importance of different physical constraints\. The terms weighted byα\\alphaandβ\\betaact locally, enforcing agreement with the reference velocity field and local physical constraints such as zero velocity inside obstacles and divergence\-free flow\. In contrast, the terms weighted byγ\\gammaandδ\\deltaoperate at a global level, enforcing periodic boundary conditions and consistency of derived macroscopic quantities, such as tortuosity\.

These parameters can therefore be interpreted as balancing local fidelity against global physical consistency\. While their values are selected through standard hyperparameter tuning, they also provide a degree of flexibility, allowing the user to emphasize specific physical properties depending on the application\.

To further investigate their influence, we performed additional experiments in which each hyperparameter was varied over approximately two orders of magnitude, while the others were kept fixed\. The tested values were:α=1,5,10\\alpha=\{1,5,10\},β=0\.1,1,10\\beta=\{0\.1,1,10\},γ=0\.01,0\.1,1\\gamma=\{0\.01,0\.1,1\}, andδ=0\.001,0\.01,0\.1\\delta=\{0\.001,0\.01,0\.1\}, where the middle values correspond to those used in the main manuscript\. Each experiment was repeated three times with different data splits, and median values are reported\.

Table[S1](https://arxiv.org/html/2605.20250#A3.T1)summarizes the results\. The table is organized in blocks corresponding to the parameter being varied, with the associated loss term highlighted for clarity\. Increasing a given hyperparameter consistently reduces the corresponding component of the loss function, as expected\. However, this does not necessarily translate into a reduction of the overall mean squared error \(MSE\), indicating non\-trivial interactions between the different loss terms\.

In many cases, adding additional loss components improves the overall prediction accuracy, which may be attributed to improved optimization and avoidance of local minima\. Notably, the parameter set used in the main manuscript \(α=5\\alpha=5,β=1\\beta=1,γ=0\.1\\gamma=0\.1,δ=0\.01\\delta=0\.01\) corresponds to one of the lowest MSE values observed\.

The interplay between local and global terms is not trivial\. For example, introducing the global periodicity term \(γ\>0\\gamma\>0\) reduces the MSE, even though it primarily enforces a global constraint\. Adding the tortuosity term \(δ\>0\\delta\>0\) slightly increases local error metrics but significantly improves the accuracy of predicted macroscopic quantities\.

It is worth noting that the same set of hyperparameters, i\.e\.α=5\\alpha=5,β=1\\beta=1,γ=0\.1\\gamma=0\.1, andδ=0\.01\\delta=0\.01, was also used in the Fourier Neural Operator experiments described in Section[S4](https://arxiv.org/html/2605.20250#A4)\. A similar trend is observed in that case: the inclusion of additional loss terms leads to a consistent improvement in all evaluated metrics\. This indicates that the beneficial effect of the proposed loss formulation is not specific to convolutional architectures, but extends to operator\-learning approaches as well\.

Overall, these results indicate that the loss function components interact in a coupled manner, and that the choice of hyperparameters reflects a trade\-off between different physical properties\. The selected parameter values represent a balanced compromise between local accuracy and global physical consistency\.

Table S1:Performance metrics for models trained with different values of loss\-function hyperparameters\. All models use the same ResNet\-101 architecture\. Median values over three independent runs with different train/validation/test splits are reported\. The table illustrates how individual loss components and prediction accuracy respond to changes in the weighting coefficients\.
## Appendix S4Fourier Neural Operators

Table S2:Comparison between the best\-performing CNN architecture \(ResNet101, see Table 2 in the manuscript\) and the Fourier Neural Operator \(FNO\) models\. Two FNO variants are considered: one trained with a simplified loss function and one trained using the full loss formulation adopted for CNNs\. In all FNO cases,α=5\\alpha=5\.In addition to convolutional neural networks, we explored the use of Fourier Neural Operators \(FNOs\) as an alternative model class\. FNOs are a recently proposed class of neural operators designed to learn mappings between function spaces, rather than between discrete imagesLi et al\. \([2021](https://arxiv.org/html/2605.20250#biba.bib3)\); Kovachki et al\. \([2023](https://arxiv.org/html/2605.20250#biba.bib4)\)\. They achieve this by performing global convolutions in the Fourier domain, enabling efficient modeling of long\-range spatial dependencies\. This makes them particularly appealing for problems governed by partial differential equations, such as fluid flow in porous media\.

In our implementation, the FNO model operates on regular grids and consists of a sequence of spectral convolution layers interleaved with pointwise transformations\. The input is augmented with spatial coordinate information, allowing the network to learn position\-dependent mappings\. The model was trained to predict the two components of the velocity field directly, using the same dataset and train–validation split as for the CNN\-based models\. After preliminary exploration of model capacity, we selected a configuration with 16 retained Fourier modes per spatial dimension, a channel width of 256, and 4 stacked spectral layers, which provided the best performance among the tested variants\.

Two versions of the FNO model were trained\. The first one employed a simplified loss function, corresponding to Eqs\. \(11–13\) in the manuscript, withα=5\\alpha=5andβ=γ=δ=0\\beta=\\gamma=\\delta=0\. The second version used the full loss formulation introduced for the CNN models, withβ=1\\beta=1,γ=0\.1\\gamma=0\.1, andδ=0\.01\\delta=0\.01\. This allows us to assess whether the additional physics\-informed terms in the loss function provide consistent benefits across different model architectures\.

The results are summarized in Table[S2](https://arxiv.org/html/2605.20250#A4.T2)\. The FNO model achieves predictive performance of the same order of magnitude as the CNN\-based approach, particularly for the velocity field reconstruction\. More importantly, the inclusion of the additional loss terms leads to a systematic improvement across all evaluated metrics\. In particular, the tortuosity prediction error decreases and the coefficient of determinationR2R^\{2\}increases when the full loss formulation is used\. A similar trend is observed for permeability\.

These observations indicate that the proposed loss function contributes to improved model performance independently of the underlying architecture\. In particular, the gains obtained by incorporating the additional loss terms are not specific to convolutional networks, but extend to operator\-learning approaches such as FNO\. This suggests that the loss formulation plays a key role in guiding the model towards physically consistent predictions\.

## References

- Fritsch and Carlson \(1980\)Fritsch, F\.N\., Carlson, R\.E\., 1980\.Monotone Piecewise Cubic Interpolation\.SIAM Journal on Numerical Analysis 17, 238–246\.[10\.1137/0717021](http://dx.doi.org/10.1137/0717021)\.
- Fritsch and Butland \(1984\)Fritsch, F\.N\., Butland, J\., 1984\.A Method for Constructing Local Monotone Piecewise Cubic Interpolants\.SIAM Journal on Scientific and Statistical Computing 5, 300–304\.[10\.1137/0905021](http://dx.doi.org/10.1137/0905021)\.
- Li et al\. \(2021\)Li, Z\., Kovachki, N\., Azizzadenesheli, K\., Liu, B\., Bhattacharya, K\., Stuart, A\., Anandkumar, A\., 2021\.Fourier Neural Operator for Parametric Partial Differential Equations\.International Conference on Learning Representations \(ICLR\)\.
- Kovachki et al\. \(2023\)Kovachki, N\., Li, Z\., Azizzadenesheli, K\., Liu, B\., Bhattacharya, K\., Stuart, A\., Anandkumar, A\., 2023\.Neural Operator: Learning Maps Between Function Spaces\.Journal of Machine Learning Research 24, 1–97\.