Simulation of near‐saturated flow in soil macropores with the lattice Boltzmann method

Direct information about the soil structure can be obtained with X‐ray computed microtomography, and the imaged macropore networks can be used as geometries in the lattice Boltzmann flow simulations. This method has not been widely applied for near‐saturated flows due to methodological issues related to diffuse‐interface two‐phase flow simulations in samples with limited resolution. Here, a simple pore‐scale lattice Boltzmann approach to simulate steady‐state water flow in partially saturated soil macropore networks that circumvents these problems is presented. The actual simulation is preceded by the determination of water–air distribution, for example, by using morphological operations. Flow through the water‐filled part of the pore space is done by using no‐slip conditions at water–solid boundaries and slip conditions at water–air interfaces. The method was tested by simulating film flow over a solid surface, and the numerical results are shown to agree with the analytical expression available for this flow geometry. The method is further tested, and its use is demonstrated with real tomographic reconstructions of clay soil samples.

depends on the water saturation level of the soil such that a decreasing saturation reduces the hydraulic conductivity (Gerke & van Genuchten, 1993).
Conventionally, hydraulic conductivity is measured in field conditions with various infiltrometry methods (Angulo-Jaramillo et al., 2000).So-called tension infiltrometers allow measurement of near-saturated hydraulic conductivity in field conditions under a small water tension by which the water saturation level can be changed in a controlled way (Ankeny et al., 1988;Clothier & White, 1981).Tension infiltrometer measurements typically use water tensions ranging from a few millimeters up to 10 cm, whereby they probe macroporosity with pore diameters exceeding ca.300 μm, which has been considered the essential pore size regime for rapid non-equilibrium flow (Jarvis, 2007).Recently, Blanchy et al. (2023) systematically collected and analyzed data on hydraulic conductivity measurements from a large number of publication sources.Their results indicate that different flow Vadose Zone Journal mechanisms are dominant at very close to saturation as compared with the dryer range.These results suggest that the largest macropores become the dominant flow mechanism in the near-saturated regime.
As measurements of hydraulic conductivity are tedious, it is necessary to develop pedotransfer functions for larger scale continuum modeling applications (Bouma, 1989).Pedotransfer functions link soil properties such as texture and organic carbon content, which are easy to measure with hydraulic conductivity (e.g., Schaap et al., 2001;Tóth et al., 2015;Zhang & Schaap, 2019).However, as near-saturated hydraulic conductivity is largely determined by the soil macropore network, they often perform quite weakly, especially for fine-textured soils with structural macropores.Pedotransfer functions for saturated and near-saturated conditions could be improved by including predictors linked to the macropore network morphology.Thus, quantitative characterization of soil macropore structure has been considered as a possibility to improve the description of hydraulic conductivity (Zhang & Schaap, 2019).
Direct information on the properties of soil structure (Anderson & Hopmans, 2013) and the soil macropore network can be obtained by three-dimensional (3D) imaging with X-ray computed microtomography (e.g., Hyväluoma et al., 2012;Larsbo et al., 2014;Schlüter et al., 2020).Soil macropore network can be segmented from imaged intact soil samples and the properties of the macropore system can be quantitively analyzed for structural characteristics such as porosity, size distribution of pores, and pore connectivity.While such simple structural properties do not directly determine the flow properties in soil pore network, the geometry and topology of pore space impose the flow friction, which results in observed hydraulic conductivity.The structural properties of soil macropore network can be associated with macroscopic flow behavior and for soil macropore systems so-called critical pore diameter has turned out to be an important parameter determining the hydraulic conductivity (Koestel et al., 2018;Schlüter et al., 2020;Soinne et al., 2023).
In addition to analyzing structural properties, macropore networks determined with X-ray microtomography can be used to obtain direct information about hydraulic conductivity by using them as simulation geometries in numerical flow simulations (Hyväluoma, Niemi, et al., 2018;Schlüter et al., 2020).One suitable approach for such simulations is the lattice Boltzmann (LB) method, which has proven to be suitable for simulating single-phase and multiphase flows in porous media (Liu et al., 2016).In the LB method, arbitrary simulation geometry can be used and typically these are obtained by 3D imaging of real porous samples, such as intact soil cores.Despite this, applications of the LB method for simulations of flow in unsaturated soils are rare as compared with singlephase saturated flows.The application of full two-phase flow simulation in porous media is restricted by the unrealistically low viscosity ratio (order of ∼10) in the most widely used

Core Ideas
• Lattice Boltzmann method was used to simulate near-saturated flow in the imaged soil macropore system.• The approach uses static water-air distribution in soil macropores that were determined by a morphological model.• Slip boundary condition at water-air interfaces was realized with a bounce-back rule from a moving boundary.
• The approach reduces the computational costs for simulating water flow in pore networks at near-saturated condition.
LB multiphase models as well as the diffuseinterface nature of these models, where the phase interface is not a sharp line but has a finite width (Liu et al., 2016).Description of the pore system of heterogeneous porous materials such as soils necessitates the use of a fairly large sample size to capture a representative volume element.Large sample size limits the discretization level that can be used in simulations due to computational restrictions (Mattila et al., 2016), and also X-ray tomography limits discretization as imaging resolution scales with sample size.Simulations of unsaturated flow in porous media are needed to gain further insight into the near-saturated flow in soil.In this note, I present an approach to studying near-saturated flow in soil with image-based LB simulations utilizing real pore geometries obtained with X-ray tomography.The proposed method uses static water distribution in soil pore space that can be determined, for example, by using an image analysis approach or a two-phase simulation of the distribution of water and air in a given porous sample.With the help of such water distribution, flow is then simulated in the water-filled part of the pore space by using the no-slip boundary condition at water-solid boundaries and the slip condition at water-air interfaces.To the best of the author's knowledge, only one model using a comparable simulation strategy has been previously presented (Zhang et al., 2016).That method is based on the so-called mirror-image method (Zhang et al., 2002) for the slip boundary treatment, while the present method relies on the bounce-back boundary condition for moving boundaries (Ladd, 1994).

Water distribution in soil macropore network
To determine the distribution of water and air in the 3D soil pore system obtained with X-ray tomography, a simple F I G U R E 1 Water-and air-filled pores in a cross-section of a soil sample.The water tension heads in the panels are 0, 1.6, 2.5, and 4.9 cm from left to right.Gray, black, and green colors denote the solid, water, and air phases, respectively.The size of each cross-section is 2.8 cm × 2.8 cm.morphology-based approach was used.All pore voxels were classified as water or air filled depending on whether the pore size related to that voxel is greater or smaller than a given threshold pore radius r.This classification of voxels was done with the morphological opening operator, that is, morphological erosion followed by morphological dilation of the set of pore voxels (Horgan, 1998).Morphological operations on pore space were done with a spherical structuring element with a radius r.The opening removes pores smaller than the threshold radius from the pore space.These pore voxels are labeled as water-filled pore space and the remainder of the pores is the air-filled pore space.Effectively, a pore voxel is considered water filled if a sphere of radius r cannot be fitted inside the pore space such that it contains that voxel.Further details of the approach used can be found in, for example, Hyväluoma, Kulju, et al. (2018).The distribution of water and air calculated for pore radius r (m) can be related to water tension h (m) and capillary pressure p c (Pa) using the Young-Laplace equation: where ρ, γ, and g are the density of water, the surface tension of water-air interface, and gravitational acceleration, respectively.In this work, the following values were used: ρ = 998 kg m −3 , γ = 0.072 N m −1 , and g = 9.81 m s −2 .Here, it is assumed that water perfectly wets the pore walls (i.e., a vanishing contact angle).Examples of the water distributions obtained with the used approach are shown in Figure 1.

Lattice Boltzmann method
The flow simulations were performed with the LB method.The basic concept, detailed theory, and application of the LB method are covered in several books and reviews (e.g., Huang et al., 2015;Liu et al., 2016;Succi, 2001;Sukop & Thorne, 2006).Here, only the aspects of the method that are necessary for the present purpose are briefly described.
LB method is based on a discrete version of the Boltzmann equation, which is solved for single-particle velocitydistribution functions.Note that all equations and quantities in this paper are written in dimensionless lattice units where the grid spacing and time step equal unity.In the LB method, fluid is modeled as fictitious fluid particles that stream and collide on a lattice, and the hydrodynamic fields (pressure, velocity) are obtained as moments of the single-particle velocity distribution function.The flow velocity obtained from the LB equation follows the Navier-Stokes equation.As the LB method operates on a cubic lattice, the pore system obtained by X-ray tomography is directly applicable to simulations, with each voxel of the image representing a lattice node in the simulation.In this work, the so-called D3Q19 model was used.This model operates on a 3D cubic lattice with 19 discrete velocities c i .In a lattice node x at time t, the state of the system is fully described by a set of particle velocity distribution functions f i (x,t), each of which gives the probability density of particles moving with velocity c i .The dynamics of the fluid is given by the LB equation: (2) Above the two-relaxation-time collision operator (Ginzburg et al., 2008) was employed, where τ e and τ o are the relaxation times for even and odd parts of the distribution function defined as (3) respectively, and direction -i is defined such that  − =   .
Here, the D3Q19 model operating in 3D cubic lattice with 19 discrete velocities was used (Qian et al., 1992).The equilibrium distribution function has the following form: where density ρ and velocity u are obtained as the zeroth and first moments of the distribution function, and   = 1∕ √ 3 is the speed of sound.Weight factors w i depend on the length of the corresponding discrete velocity c i (Qian et al., 1992).Even and odd equilibrium distributions  ()  and  ()  are defined as in Equation (3).The kinematic viscosity ν is determined by the even relaxation time ν =  2  (  − 1∕2) and the even and odd relaxation times are related by the "magic formula" (d'Humiéres & Ginzburg, 2009): In this work, the LB method was used to simulate flow in the water-filled part of the soil macropore system.At water-solid boundaries, the no-slip boundary condition was enforced with the bounce-back scheme, where distribution functions colliding with a solid node are bounced back to the fluid domain, that is,   =  − .The handling of slip boundary conditions at the water-air interface is described in Section 2.3.
The flow was driven by a body force in the z direction (Guo et al., 2002).After obtaining the steady flow field, the permeability k of the sample was determined using the Darcy law: Here, q z is the volumetric water flux in the flow direction, μ = ρν is the dynamic viscosity, and ∂p/∂z is the pressure gradient in the flow direction.The hydraulic conductivity can then be calculated from permeability as K = k/(ρg), where ρ is the density of water and g is the gravitational acceleration.The simulations here were performed in laminar flow regime where Darcy's law is valid.

Boundary conditions at the water-air interface
At water-air interfaces, the bounce-back boundary condition was modified.As the viscosity ratio between water and air is large and the shear stress at the interface is negligible, a no-slip boundary condition can be used.Here, the no-slip condition is realized by using the bounce-back from a moving boundary (Ladd, 1994): Here, u B is the local velocity of the water-air interface, x is a water node adjacent to the interface, and velocity vector c -i points from the water domain toward the air-filled part of the pore space.The velocity of the boundary is approximated by the water velocity in node x.However, since in this approach the water-air interfaces are rigid, the velocity component normal to the water-air interface is subtracted from the boundary velocity, that is, Here, n is the surface normal vector for the water-air interface, which is calculated by using Sobel operators as described by Hyväluoma, Niemi, et al. (2018).Without subtraction of the normal velocity component, simulations driven by a forcing term could easily develop numerical instabilities.

NUMERICAL EXPERIMENTS
As a first benchmark, the above-described approach was tested by simulating the flow of a water film over a solid wall (Figure 2a).At the water-solid boundary, a no-slip boundary condition was enforced, whereas at the water-air interface, slip boundary condition was applied.The velocity field in this simple flow system can be obtained analytically by solving the Navier-Stokes equations with boundary conditions u x = 0 at the solid wall and   ∕z = 0 at the air interface.The solution reads as where D is the film thickness, and coordinates x and z are as defined in Figure 2a.Simulations were performed by varying the film depth from 6 to 100 lattice spacings, and the pressure gradient was chosen so that the flow remained in the low-Mach-number regime (i.e., velocity was clearly lower than speed of sound).An example velocity profile from simulation and the corresponding analytical velocity profile are shown in Figure 2b.Simulation error was assessed by the L 2 error norm, defined as where u x and U x are the simulated and analytical flow velocity, respectively.In Figure 2c, shown is the L 2 error norm calculated for different film thicknesses, which indicates that the error converges as film thickness (i.e., discretization level) increases.
As a second example, flow through pore systems in real clay soil samples imaged by X-ray tomography was simulated.The sample size in these simulations was 941 × 941 × 667 voxels (2.8 × 2.8 × 2.0 cm 3 ), and the isotropic voxel size was 30 μm.Three images with different X-ray visible macroporosities were selected from a larger set of samples collected and imaged in a previous project, which are hereafter denoted as samples S1, S2, and S3.Image analysis results for the full set of samples have been reported elsewhere, along with detailed information about imaging and image processing (Soinne et al., 2023).Segmented images, with every voxel classified either as a pore or solid voxel, were used in the present simulations.For each sample, water distribution within the pore space was determined for seven threshold pore radii ranging from 60 to 900 μm, which corresponds to a water tension head between 24.5 and 1.6 cm.These are also typical tensions used in tension infiltrometer measurements (e.g., Carey et al., 2007;Jarvis & Messing, 1993).Examples of the water distributions for sample S2 for different tensions are shown in Figure 1.The X-ray visible macroporosities of the samples were 3.5%, 7.8%, and 10.4% for samples S1, S2, and S3, respectively.The water-filled volume fraction of the pore space for each pressure head is shown in Figure 3a, where the differences between samples, especially close to saturation, indicate varying pore size distributions.For each sample and water saturation level, the flow was driven by a body force representing the gravity in the z direction until the steady-state flow was approached.The Darcy law was then applied to obtain hydraulic conductivity for each sample and water tension.Results are shown in Figure 3b as the relative hydraulic conductivity using the saturated hydraulic conductivity as reference value,   (ℎ) = (ℎ)∕  .As for water-filled porosity, the relative permeabilities also differ, especially in the wet-end part of the curve (h < 5 cm).The slight increase in the relative hydraulic conductivity observed between the two highest tensions is probably due to insufficient discretization of the smallest pores.

DISCUSSION
An LB method for simulating flow through a partially saturated complex pore structure was presented.This model is very simple and easy to implement, and the generic formulation of the boundary conditions at water-air interfaces by a bounce-back rule from moving wall ensures that all possible geometric configurations can be handled.Therefore, it can be easily added to existing LB implementations.The results of the benchmark case with film flow over a solid surface agreed with the analytical solution available for this case.However, the error showed only first-order convergence.This is weaker than the formal second-order convergence of the method, which is nevertheless dependent on the boundary treatment (Chen et al., 2020;Holdych et al., 2004).This difference likely stems from the zeroth-order extrapolation of the water-air interface velocity from the neighboring water node.Previously, a comparable reduction in error convergence has been observed in LB simulations of wall shear stress at solid surfaces (Hyväluoma, Kulju, et al., 2018).In that case, the use of higher order extrapolation schemes has been shown to offer possibilities for improvement (Kang & Dung, 2014).However, in the present case, such approaches would complicate the use of the method in arbitrary complex geometries.As a more complex test case, flow through pore networks obtained with X-ray tomography was simulated.While the number of samples used in this demonstration was small (three samples), the results emphasize that the influence of macroporosity on near-saturated hydraulic conductive requires careful consideration.Blanchy et al. (2023) observed in their meta-analysis that near-saturated hydraulic conductivities measured at tensions between 4 and 10 cm were mutually relatively well correlated.In contrast, correlations between saturated hydraulic conductivity and near-saturated hydraulic conductivities at tensions smaller than 4 cm were only poorly correlated to each other.The present results show a similar tendency.These observations emphasize the importance of considering the near-saturated regime in addition to saturated hydraulic conductivity when pore-scale simulations with Xray images are performed.Saturated hydraulic conductivity has been found to relate to the critical pore diameter (i.e., the diameter of the largest sphere that can pass through the pore system in the vertical direction) via a power law dependence (Koestel et al., 2018;Soinne et al., 2023).This emphasizes the role of the largest continuous flow paths in the sample, whereas the limited sample height in X-ray-based studies can lead to unrealistically high values of saturated hydraulic conductivity.In line with these considerations, it has been observed that the saturated hydraulic conductivity is not a good calibration point for the Mualem-van Genuchten model because it is only weakly correlated to the near-saturated hydraulic conductivities (Schaap & Leij, 2000;Schaap & van Genuchten, 2006;Vogel et al., 2001).This is also related to the non-convergence of the integrals in the conductivity functions for realistic values of the van Genuchten parameters (Durner, 1994;Ippisch et al., 2006;Madi et al., 2018).
While the water-air distribution in the 3D pore system imaged with X-ray tomography was determined with the fairly simple morphological approach, other approaches can be selected as well.The selection may depend on which physical process has been used to study the water distribution or what are the requirements for the computational cost.One possibility would be to use two-phase LB simulations to determine the water distribution by simulating phase separation in the imaged pore system (e.g., Genty & Pot, 2013;Pot et al., 2015).However, as such an approach would require heavy computations in addition to the actual flow simulation, it would have been impractical for the present purpose but can be useful if, for example, a non-zero contact angle is considered.Another option could be to utilize more sophisticated mathematical morphology-based method, such as the method developed by Hilpert and Miller (2001) for simulating drainage of a porous system.The present approach is not applicable for dynamic features such as intermittent flows (Tofteng et al., 2002) or airentrapment processes (Snehota et al., 2015), but consideration of such phenomena requires using a full two-phase simulation model.

CONCLUSION
The proposed simulation approach provides a tool to analyze the effects of soil macropore network architecture on nearsaturated hydraulic conductivity.The method presented is a simple extension to the LB method, which allows simulations of near-saturated hydraulic conductivity but is limited to static air-water distributions.Therefore, the method is not capable of handling flow processes in soil with dynamic air-water interfaces.Such flow modes require computationally more intensive full multiphase simulations.However, the method can be a useful tool to relate structural properties obtained by X-ray tomography with flow properties.The method can be particularly applicable in studies where nearsaturated hydraulic conductivity needs to be determined for a larger number of X-ray tomography images, as the computational cost is similar to single-phase flow simulations and the method avoids problems related to diffuse-interface methods.While this work is focused on macropores, which are a typical pore size regime for X-ray tomography studies, the method is suitable for other pore size regimes and sample sizes as well.

AU T H O R C O N T R I B U T I O N S
Jari Hyväluoma: Conceptualization; funding acquisition; investigation; methodology; writing-original draft.

A C K N O W L E D G M E N T S
This work was funded by The Finnish Research Impact Foundation.CSC-IT Center for Science is acknowledged for computing time through project #2005253.

C O N F L I C T O F I N T E R E S T S T A T E M E N T
The author declares no conflicts of interest.

F
Benchmark test with water film flowing over a solid surface.(a) Schematic picture of the test geometry.(b) The simulated velocity profile for the film thickness of 30 lattice spacings is shown with markers, and the analytical velocity profile is shown with solid line.(c) L 2 error norm as a function of film thickness is shown with markers.The solid line shows first-order error convergence (L 2 ∝D −1 ) as a guide for the eye.Dimensionless lattice units (lu) are used.

F
I G U R E 3 (a) Water-filled porosity of the three imaged soil samples for varying water tension.(b) Relative hydraulic conductivity determined from flow simulations for samples with varying water content.Relative hydraulic conductivity is calculated by using the saturated hydraulic conductivity as reference value.