Evolution of Water Retention Characteristics in Bio-Geochemically Altered Unsaturated Soils: A Pore-Scale Study

. Biogeochemical processes in subsurface can dramatically alter the behavior of multiphase fluid flow and the hydrodynamics of porous media by bio-clogging. There exist analytical solutions such as soil suction – saturation equations which can be used to predict the water retention curve relations in unsaturated soils. However, due to the complexity of various biogeochemical products, their pore-scale behavior and their interplay with pore structure, such analytical solutions would not provide accurate predictions. In this study, a large database of pore-networks were generated by adjusting the statistical and spatial pore and tube size distribution of the networks resembling various levels of bio-clogging. Numerical simulations including the evolution of pore structure, water retention relationship, and air invasion dynamics during desaturation were explored. The numerical simulations verified that local pore-clogging leads to the development of isolated pore clusters and impermeable zones. The evolution of impermeable zones results in the formation of preferential flow paths towards the mobile zones for the multiphase flow problems. Using parallel computation, the critical predictors of water retention curves in bio-clogged porous media are found.


Introduction
The fundamental expressions for unsaturated soil behavior are anchored in capillarity and water saturation relation Pc-Sw, typically known as a water retention curve (WRC). The water retention curve essentially captures pore-scale characteristics of the porous media and the gas-fluid-mineral interactions. This curve can be used as an indicator to hydraulic conductivity [1], soilwater storage [2], and soil stiffness, strength, and volume changes [3][4][5][6]. The relative fluxes of liquid and gaseous phases through porous media are usually captured by relative water and gas permeability krw and krg. Both the water retention curve and relative permeability are the most critical characteristics to characterize multiphase flow in porous media.
Biogeochemical processes in subsurface including biomineral, biofilm, and biogas formation can dramatically alter the behavior of multiphase fluid flow and the hydrodynamics of porous media. However, acquiring water retention curve for bio-geochemically altered soils in the laboratory is extremely timeconsuming and also involves experimental limitation to dictate various bio-clogging mechanisms in soils.
Widely used models of water retention curve and relative permeability [7][8][9][10][11][12] contain multiple fitting parameters. The appropriate selection of these * Corresponding author: nmahabadi@uakron.edu parameter values is critical to enhance the prediction of unsaturated behavior of bio-clogged soils.
In this study and as a first step towards a clear understanding on the evolution of porous media during bio-geochemical processes and its impacts on WRC, the multiphase flow of gas invasion into water saturated pore space of bio-clogged soils are obtained using Pore-Network simulation. A large database of pore-networks were generated by adjusting the statistical and spatial pore and tube size distribution of the networks resembling various levels of bio-clogging. Numerical simulations including the evolution of pore structure, water retention relationship, and air invasion dynamics during desaturation were explored. The numerical simulations verified that local pore-clogging leads to the development of isolated pore clusters and impermeable zones. The evolution of impermeable zones results in the formation of preferential flow paths towards the mobile zones for the multiphase flow problems. Using parallel computation, the critical predictors of water retention curves in bio-clogged porous media are found.

Pore-Network Model
In order to find the effect of pore-scale characteristics of soi medium during bio-clogging and its impacts on water retention curve, a significant series of computational fluid dynamic (CFD) simulations is required. However, due to the complexity involved in natural porous systems, direct numerical methods such as Finite element Methods which requires mesh generation and solving partial differential equations on the entire domain geometry are prohibitively time consuming and expensive. To overcome this limitation, the numerical simulations performed in this study are formed around several numerical algorithms developed by Mahabadi et al. 2016 [13], Mahabadi et al. 2020 [14] which aims to study complex coupled processes in highly heterogenous porous media of natural soil and rock systems. The basis of these algorithms is built on the concept of Pore Network Models (PNM). A PNM consists of interconnected networks of nodes (Pores) and tubes (throats) which can be either directly extracted from natural sediments or numerically generated to simplify and satisfy the geometry of the original natural porous network. Using the PNM, one can efficiently simulate the flow of a viscous fluid in the entire network by considering the conservation of mass and solving Hagen-Poiseuille equation which can significantly reduce the computational cost of numerical simulations. The other advantage is the ability of pore network concept to simulate a large database of different porous media with varying pore-scale characteristics by tuning the statistical information.

Pore-Network Generation
A large number of pore networks (over 100,000) are generated by adjusting statistical pore-scale characteristics of networks to simulate the multiphase flow during gas invasion and calculate the WRC parameters for a wide variety of bio-clogged porous media.
The baseline of pore network models is generated using a 25 by 25 pore lattice system, with a poreconnectivity (Coordination Number, CN) of 8, in which each pore is connected to all its neighboring pores (CN=8). The mean pore and tube (throat) size and their size distribution (standard deviation of pores and tubes, σpore and σtube) are varied to cover a large range of pore size distributions. Particularly, as the capillarity is driven by the movement of water-gas meniscus in pore throats, the mean size of tubes (Rtube) is randomly selected between 1 nm to 1 cm resembling a wide range of soil media from clayey to gravel soils. The standard deviation of pore and tube size distributions (σ) are also varied to generate porous systems ranging from highly uniform networks (σ=0) to highly heterogenous porous media (σ=1).

Clogging in PNMs
In this study, clogging of porous media is manipulated by random elimination of tubes from the pore networks. Initially, the connectivity number of the networks was 8 (CN=8), then networks are subjected to randomly elimination of tubes to simulate clogging. The clogging of each pore network is generated from 0% and continued until fluid flow terminates in the pore networks (in some cases the final clogging is up to 80%). As the clogging in networks cause an evolution on the pore scale characteristics of the database, after each elimination, the statistical information of clogged networks such as the average connectivity number, impacted pore and tube size, etc are calculated. For example, the clogging in the networks caused the average CN to range from CN=8 (baseline non-clogged network) to less than CN<2 (for highly clogged pore networks).

Gas Invasion (Desaturation)
In the second stage of simulation, it was assumed that the network is fully water saturated. The capillary pressure was calculated for each tube based on Young-Laplace equation: Where Pc is the capillary pressure in the tube that has to be replaced by gas, Rtube is the radius of tube, γ is the surface interfacial tension (0.072 N/m), θ is the wetting angle of the liquid on the surface of the capillary, which is 0° in this study to represent a completely water-wet surface [15]. Based on the equation (1), capillary pressure for each tube was calculated. The process of desaturation were simulated Gas invasion is enforced at tubes located on the inlet boundary by gradually increasing the gas pressure.
Quasi-static gas invasion controlled by capillarity is assumed so the viscosity effect can be disregarded. For each invasion step, the candidate tubes for gas invasion (the water-filled tubes in the connection to gas phase) are found where the capillary pressure is calculated for each candidate. If the gas pressure is higher than the capillary pressure of the water-filled tube candidate, the gas can invade the tube. In this study, we assumed that the water pressure inside of the tubes is negligible (Pw=0). With pressure increases, gas starts invading through the largest tube on the inlet boundary. And as capillary pressure increases, gas invades more waterfilled tubes, leading to a decrease in water saturation; thus, the water retention curve is computed until water drainage stops.

Gas Invasion Dynamics
Here, we explain the impact of clogging on hydrodynamics of the two-phase flow when air invades into the fully saturated pore networks. Figure 2 and 3 shows the development of air invasion for a uniform (σ=0) and non-uniform (σ=0.2) pore network respectively. As shown in figure 2, in a non-clogged (CN=8) uniform pore-network, there is a clear homogeneous front of the gas phase which displaces the water in the porous media. In this case, all the pores filled with water are displaced with air where there is no residual wate saturation at the end of the invasion process. However, as the clogging increases (CN decreases), air invasion patterns exhibit higher degrees of heterogeneity following by higher residual water saturation in the system. Although figure 2 shows the air invasion for a completely uniform porous media (σ=0), it is clear that higher percentage of clogging will transform the hydrodynamics of air invasion by formation of preferential pathways of the air phase. In this case, there is significant volume of water pores cannot be displaced by air. Fig. 2. Air invasion in a uniform pore network (σ=0) generated with various connectivity numbers (from CN=8 for nonclogged to CN=2 for highly clogged porous medium). Figure 3 shows the invasion process for a non-uniform network (σ=0.2). In this case, clogging is not the only feature that defines the mechanism of air invasion patterns. The heterogeneity of porous media also plays a critical role to explain the behavior of air invasion into the saturation system. Fig. 3. Air invasion in a non-uniform pore network (σ=0.2) generated with various connectivity numbers (from CN=8 for non-clogged to CN=2 for highly clogged porous medium).

Water Retention Curves
The water retention curves are generated by calculating the saturation at each gas invasion iteration associated with each gas pressure required to overcome the capillary pressure of the candidate tubes in the network. Several analytical models are suggested to predict water retention curves. van Genuchten (1980) model has been used widely due to its simplicity in which the capillary pressure-water saturation (Pc-Sw) relation can be generated using three variables of air entry pressure (P0), residual water saturation (Sw) and the fitting parameter of m.
Lower m-value manipulate steeper Pc-Sw curve ( Figure  3), resembling a wider pore size distribution of sediment. Figure 4 shows typical water retention curves with varying fitting parameters m, using the van Genuchten model.   The data presented is about 3000 simulations (out of 100,000) for the pore networks with the air entry value (P0) ranges between 10 kPa to 100 kPa. Different colors contribute to different m-fitting values (slope of the curves).
Using MATLAB curve fitting tool, the simulation results for the entire 100,000 pore networks (100,000 water retention curves) were fitted via van Genuchten equation (equation 2) to find the three input parameters of P0, Srw and m.

Stochastic analysis
The entire database of 100,000 simulation of water retention curves is analysed to explore the impact of clogging on the behavior of multiphase flow hydrodynamics during air invasion. Particularly, we searche for the most important variables which impact the behavior of WRC. To do this analysis, we use "Permutation Importance" analysis. Permutation feature importance measures the increase in the prediction error of the model after we permuted the feature's values, which breaks the relationship between the feature and the true outcome. We measure the importance of a feature (e.g. pore connectivity, pore space heterogeneity, clogging, etc) by calculating the increase in the model's prediction (e.g. WRC behavior) error after permuting the feature. A feature is "important" if shuffling its values increases the model error, because in this case the model relied on the feature for the prediction. A feature is "unimportant" if shuffling its values leaves the model error unchanged, because in this case the model ignored the feature for the prediction. The permutation feature importance measurement was introduced by Breiman (2001)43 for random forests.
The following sections illustrates the impact of important features on the behavior of WRC. Figure 6 shows the impact of clogging on the evolution of porous media in terms of pore connectivity number (CN) and percentage of the isolated pores. It is found that there is a linear relationship between pore connectivity and clogging percentage (Figure 6-a). In this case the 3 rd important feature is found to be the isolated pores (colorbar). Now, we plot the relationship between these three features (CN, clogging and isolated pores) in a different way (Figure 6-b). In this case, the Higher the clogging, lower the pore connectivity number, and higher the probability of more isolated pores in the system. However, it is interesting that the relationship between the pore connectivity number and isolated pores is in the form of an exponential function. Particularly, a significant growth of isolated pores (%) is found when CN<3. In this case (CN<3), the formation of preferential air flow pathways (fringes) cause a large volume of pores being immobilized.

Impact of clogging on Residual Water Saturation
Now, we investigate the impact of various features on the WRC. For residual water saturation (Srw), the most important features found to be CN and clogging (%). Figure 7 shows the relationships between these variables. Again, the impact of clogging is significant on the formation of isolated pores which in turn increases the residual water saturation. Higher the clogging, lower the pore connectivity number, and higher the residual water saturation. Here, the variation of Srw for a given CN cannot be explained by only the clogging (5). It is found the heterogeneity of the porous media (σ) cause this variation.  . 7. Impact of clogging on the pore connectivity number and the residual water saturation (Srw)

Impact of clogging on fitting parameter m
As explain in figure 4 and 5, different porous media would not only exhibit various air entry value (P0) and residual water saturation (Srw), but also different Pc-Sw slope. The slope of Pc-Sw can be explained by the fitting parameter of m in van Genuchten model. Using the permutation importance analysis of the forest model, the important features for fitting parameter of m are found to be the clogging (%) and the heterogeneity of the porous media (σ). Figure 8 shows the relationship between these three variables. For highly uniform networks (σ<0.1), the m-value is close to 1 (resembling a flat WRC curve) regardless of the change in clogging.
In this case, even with high percentage of clogging, there is a negligible decrease in the m-value. However, as the heterogeneity of pore space increases, there is a higher probability of variation in the m-value. In general, higher the heterogeneity (higher σ), lower the m-value. This is can be supported by the fact that higher σ, generates higher variation in the pore and throat size distribution, which in turn requires a large diversity of different air pressures to allow invasion into the pore networks. Fig. 8. Impact of clogging on the van Genuchten fitting parameter of m and heterogeneity of the porous media (σ).

Impact of clogging on air entry value
The key features affecting the air entry value (P0) is found to be size of pores, and the heterogeneity of porous media. Particularly, the larges pore throats of each network (Rmax) describe the air entry value of that network. Figure 9 shows the relation between Rmax and P0. By incorporating equation 1, we also calculate the theoretical air entry value using the Rmax for each network. As shown in figure 9, the lower boundary of the database completely matches with this theoretical line. However, as the heterogeneity of the system increases (higher σ), there will be a higher deviation from the theoretical line. This is due to the fact that the largest pore throat of each network is not necessarily the one associated with the air entry value (first throat invaded by air). If the network is completely uniform, then all the throats have same size. Therefore, Rmax defines P0. However, if for higher σ, the throat size distribution is higher where the largest throat may be located somewhere in the middle of the network (cannot be invaded early in the invasion process).

Fig. 9.
Relationships between air entry value (P0), the maximum size of throats (Rmax) and heterogeneity of the porous media (σ).

Conclusion
Biogeochemical processes in subsurface can dramatically alter the behavior of multiphase fluid flow and the hydrodynamics of porous media by bioclogging. There exist analytical solutions such as soil suction -saturation equations which can be used to predict the water retention curve relations in unsaturated soils. However, due to the complexity of various biogeochemical products, their pore-scale behavior and their interplay with pore structure, such analytical solutions would not provide accurate predictions. In this study, a large database of pore-networks were generated by adjusting the statistical and spatial pore and tube size distribution of the networks resembling various levels of bio-clogging. Numerical simulations including the evolution of pore structure, water retention relationship, and air invasion dynamics during desaturation were explored. The numerical simulations verified that poreclogging leads to the development of isolated pore clusters and impermeable zones. The evolution of impermeable zones results in the formation of preferential flow paths towards the mobile zones for the multiphase flow problems. Using parallel computation, the critical predictors of water retention curves in bioclogged porous media are found.