Statistical Effective Diffusivity Estimation in Porous Media Using an Integrated On-site Imaging Workflow for Synchrotron Users

Transport in porous media plays an essential role for many physical, engineering, biological and environmental processes. Novel synchrotron imaging techniques and image-based models have enabled more robust quantification of geometric structures that influence transport through the pore space. However, image-based modelling is computationally expensive, and end users often require, while conducting imaging campaign, fast and agile bulk-scale effective parameter estimates that account for the pore-scale details. In this manuscript we enhance a pre-existing image-based model solver known as OpenImpala to estimate bulk-scale effective transport parameters. In particular, the boundary conditions and equations in OpenImpala were modified in order to estimate the effective diffusivity in an imaged system/geometry via a formal multi-scale homogenisation expansion. Estimates of effective pore space diffusivity were generated for a range of elementary volume sizes to estimate when the effective diffusivity values begin to converge to a single value. Results from OpenImpala were validated against a commercial finite element method package COMSOL Multiphysics (abbreviated as COMSOL). Results showed that the effective diffusivity values determined with OpenImpala were similar to those estimated by COMSOL. Tests on larger domains comparing a full image-based model to a homogenised (geometrically uniform) domain that used the effective diffusivity parameters showed differences below 2 % error, thus verifying the accuracy of the effective diffusivity estimates. Finally, we compared OpenImpala’s parallel computing speeds to COMSOL. OpenImpala consistently ran simulations within fractions of minutes, which was two orders of magnitude faster than COMSOL providing identical supercomputing specifications. In conclusion, we demonstrated OpenImpala’s utility as part of an on-site tomography processing pipeline allowing for fast and agile assessment of porous media processes and to guide imaging campaigns while they are happening at synchrotron beamlines. Supplementary Information The online version contains supplementary material available at 10.1007/s11242-023-01993-7.

1 Effective Parameter Estimation 1.1 Mathematical Homogenisation Method -full derivation Formal effective parameter estimation for different physical phenomena in porous media has been well established [1].This study focuses on obtaining effective estimates for diffusivity initially considering an image based domain.
A key assumption that will be followed in this work is that there exist a separation of scales associated with the small and large scale.The small scale presents all of the imaged detail in the pore space and is assumed to be periodic in the different directions.The large scale can take the estimated effective (or homogenised) diffusivity and generate estimates that are sufficient for describing the large transport scale behaviour without explicitly considering all of the pore scale detail.
Considering the diffusion equation: where c [mol m −3 ] is a concentration of a transporting solute, t [s] is time, D [m 2 s −1 ] is the diffusivity, x is the spatial coordinate associated with regions in the domain volume Ω (Figure 1 A.), and ∇ is the spatial differentiation operator.Assuming the dimensional scaling of the different terms is where L X [m] is the large length scale, τ [s] is the time scale, and c is the average initial concentration of solute.The non-dimensional equation becomes: We assume that the domain Ω = n i Ω i , where each Ω i is a periodically repeating sub-domain (or cell) (Figure 1 A-B.).As such, the boundary condition for periodicity is given as: where ∂Ω i is the outer boundary of a cell (Figure 1 B.).We assume that solutes cannot move across internal boundaries between pore-space and solids: where Γ i represents the internal boundaries between solids and pore-space (Figure 1 B.).

Multi-scale expansion
The dimensional differential operator can be expanded into a large scale estimate and a small scale corrector: where X and ξ are the large and small scale spatial variable respectively.The dimensional scaling of the differential operators is given by: Considering the relationship between the large scale and small scale as L X >> l ξ , we say that the ratio between the scales is given by Lx = where is a small number.Our dimensionless differential operator can be expressed as: We then assume that the concentration can be expanded as: Substituting in eqs.8 and 9 into eq.3, we obtain the following equation on the cell: with the associated no flux conditions: and We note that given the periodicity of the cell, it suffices to consider D = D(ξ).
Considering only the terms in −2 , we have the following system of equations: We consider the domain integral of c (0) (∇ ξ •(D∇ ξ c (0) )) and use Green's identity to expand the integral and obtain: Ωi Due to the domain equations and boundary conditions (eq 13), only the last term remains.Thus, eq. 14 becomes: This can only be true if either D = 0 or ∇ ξ c (0) 2 = 0. Assuming the non-trivial case (i.e.D = 0), we can deduce from this that c (0) (t, ξ, X) = c (0) (t, X), thus not varying on the small scale.

Order O( −1 )
Considering the terms in −1 , we have the following system of equations: From eq 15, we can simplify the equation to: We represent the gradient term as: where êi are the unit vectors in the X i directions, and we expand eq.17 as: We now consider the following system of equations who's solutions will act as small scale correctors to our homogenised system: where χ i is the small scale corrector in the i th direction.If we multiply both sides of the eq.20 by ∂Xi , we obtain: Note that by eq. 15, ∂Xi is invariant to differentiation by the ∇ ξ (•) operator.
As such, if we take the summation of eq.21 for the different i th directions, we obtain: By substituting eq.22 into eq.19, we can see the relationship between the two scales manifests as: where c(1) vanishes under differentiation by the ∇ ξ (•) operator.

Order O( 0 )
We consider lastly the terms on the order O( 0 ): We consider the domain equation inside of the integral.As such, we have: By divergence theorem, we can see that the first term becomes: and by the boundary conditions in eq.24, we can see: Thus, eq. 25 can be re-written as: By substituting in eq.23 into eq.28, we get: where I is the n sized identity matrix and Given that c (0) is does not vary on Ω i and ∇ X (•) is separate from ξ, eq.29 can be re-written as: were Ω i = Ωi dξ.Thus, the large scale homogenised diffusion equations can be represented as: where is the effective diffusivity parameter.

1D analytic solution used for validation
To validate the results of our effective diffusivity parameter, we ran numerical simulations with the full heterogeneous geometry duplicated in a horizontal array and compared the results to a 1 D analytic solution of the same length scale uniform geometry considering the effective diffusivity.The 1D model considered the following equations: We define a variable ĉ(0) = c (0) − c ∞ , where we can re-write our equation as: We define ĉ(0) (t, X) = Y (X)T (t), which allows us to use separation of variables.We can re-write the domain equation of eq 34 as: which can then be expressed as: where λ is a constant.This generates two ordinary differential equations, one and one in time: Considering the Boundary value problem in eq.37, our general solution will take the form: Given the first boundary condition in eq.37, we know that A = 0.If we differentiate the solution in eq.39 and use the second boundary condition in eq.37, we get the expression: which implies Therefore, we have Returning to eq. 38, we know that the general solution will be of the form: Our solution can be a super position of solutions, thus ĉ(0) (t, X) = Given the initial condition of eq.34, we know: Thus, c n = − 4c∞ π(2n+1) .Since our full equation is c (0) = c ∞ + ĉ(0) , we have our analytic homogenised estimate as: e −D ef f ( π 2L (2n+1)) 2 t sin ( π 2L (2n + 1)X)). (47) 3 Influence of surface reactions on homogenisation While the scope of this study was to investigate the influence of physical geometric impedance on mass transport, we know that the homogenisation

Fig. 1
Fig.1Generating effective parameter estimate based on XCT images.A -Representation of a large scale pore space Ω comprised of periodic subdomains Ω i .The length scale of the large domain is represented by L X while the length scale of the subdomain is represented by l ξ .B -The full image based geometry representing the pore space within a porous media (soil pore space in the given example).C -Determining the representative elementary volume sufficient for obtaining reasonable effective parameter (i.e.diffusivity) estimates.D -Generating a region of interest for a given size at a random location in the modelled domain.