Next Article in Journal
Distribution Characteristics and Source Analysis of Microplastics in Urban Freshwater Lakes: A Case Study in Songshan Lake of Dongguan, China
Next Article in Special Issue
Agricultural Productive Carrying Capacity Improve and Water Optimal Allocation under Uncertainty Based on Remote Sensing Data in Lancang County, Southwest China
Previous Article in Journal
Study of Jingjiang Beach Morphodynamics in the Tidal Reach of the Yangtze River
Previous Article in Special Issue
Effect of Different Water Treatments in Soil-Plant-Atmosphere Continuum Based on Intelligent Weighing Systems
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

How Heterogeneous Pore Scale Distributions of Wettability Affect Infiltration into Porous Media

1
iES Landau, Institute for Environmental Sciences, Geophysics, University of Koblenz-Landau, Fortstraße 7, 76829 Landau, Germany
2
Institute for Bio- and Geosciences, IBG-3 Agrosphere, Forschungszentrum Jülich GmbH, Wilhelm-Johnen-Straße, 52428 Jülich, Germany
3
Institute of Building Materials and Concrete Structures (IMB), Karlsruhe Institute of Technology (KIT), 76131 Karlsruhe, Germany
4
Materials Testing and Research Institute Karlsruhe (MPA), Karlsruhe Institute of Technology (KIT), 76131 Karlsruhe, Germany
5
Physics of Soils and Terrestrial Ecosystems, Department of Environmental Systems Science, ETH Zürich, Universitätstraße 16, 8092 Zürich, Switzerland
6
Modelling and Numerics, Department of Mathematics, University of Erlangen-Nürnberg, Cauerstraße 11, 91058 Erlangen, Germany
7
Institute of Crop Science and Resource Conservation (INRES), University Bonn, Karlrobert-Kreiten-Straße 13, 53115 Bonn, Germany
*
Author to whom correspondence should be addressed.
Water 2022, 14(7), 1110; https://doi.org/10.3390/w14071110
Submission received: 24 February 2022 / Revised: 20 March 2022 / Accepted: 28 March 2022 / Published: 30 March 2022
(This article belongs to the Special Issue Numerical Modeling of Water Flow, Nutrients and Sediment Transport)

Abstract

:
Wettability is an important parameter that significantly determines hydrology in porous media, and it especially controls the flow of water across the rhizosphere—the soil-plant interface. However, the influence of spatially heterogeneous distributions on the soil particles surfaces is scarcely known. Therefore, this study investigates the influence of spatially heterogeneous wettability distributions on infiltration into porous media. For this purpose, we utilize a two-phase flow model based on Lattice-Boltzmann to numerically simulate the infiltration in porous media with a simplified geometry and for various selected heterogeneous wettability coatings. Additionally, we simulated the rewetting of the dry rhizosphere of a sandy soil where dry hydrophobic mucilage depositions on the particle surface are represented via a locally increased contact angle. In particular, we can show that hydraulic dynamics and water repellency are determined by the specific location of wettability patterns within the pore space. When present at certain locations, tiny hydrophobic depositions can cause water repellency in an otherwise well-wettable soil. In this case, averaged, effective contact angle parameterizations such as the Cassie equation are unsuitable. At critical conditions, when the rhizosphere limits root water uptake, consideration of the specific microscale locations of exudate depositions may improve models of root water uptake.

1. Introduction

Despite several decades of intensive research in the field of soil physics, the reliable prediction of water flows in soil remains a major challenge. This is mainly due to the scale and time-dependent heterogeneity characteristic of soils, which can be of a structural, physical, chemical, and biological nature [1,2]. Especially soil organic matter (SOM) plays an important role in this context. Due to its numerous input possibilities as well as highly varying degrees of degradation and mineralization, SOM can be heterogeneously distributed across scales. Their importance for water transport in soil is partly due to their strong influence on the wettability of soil surfaces [3]. The contact angle ( CA ), which is present at the three-phase boundary between liquid, gas and soil surface, serves as a characteristic quantity for the wettability [4]. A CA of zero is referred to as wettable (hydrophilic), in the range 0 ° < CA < 90 ° as reduced wettability and CA 90 ° as non-wettable (hydrophobic). In soil, hydrophilicity is mainly caused by mineral compounds [5], while hydrophobicity is mainly caused by non-polar organic compounds in the form of interstitial particulate material and coatings of mineral surfaces [6,7]. It should be noted that most soils are neither completely hydrophilic nor hydrophobic, but exhibit reduced wettability [8,9,10]. Within soil science, wettability or soil water repellency (SWR) is an important factor influencing soil erosion, surface and subsurface hydrology, and plant growth [3]. The significance of SWR for soil hydrology is based on the fact that it promotes preferential flow, surface runoff and the generation of heterogeneous moisture patterns, and thus above all the hydraulic heterogeneity of soil [11]. Despite the impact of SWR on soil hydraulics, still little is known about its spatial variability on different scales [11]. To quantify the effects of SWR on soil hydraulics, information on the specific microscopic distributions may be require. The spatial distributions depend on the specific processes which create the SWR like metabolic products of microorganisms, root and fungal exudates and their degradation products [11,12].
The rhizosphere is a narrow layer around the roots, which is a hotspot for microbial activity dominated by interactions between plant roots and soil that cause heterogenous spatial and temporal distributions of soil organic matter leading to varying SWR [13]. The pore-scale influence of SWR distributions on water transport in rhizosphere has not yet been fully investigated. Due to root growth and the release of plant root exudates, rhizosphere differs significantly in physical, chemical and biological properties from the surrounding bulk soil [14,15,16,17]. One of the most pertinent root exudates is mucilage, a hydrogel which can significantly alter the hydraulic properties of the rhizosphere [18,19,20]. It is characterized by a large water holding capacity of up to 1000 times of its own weight depending on the type [21,22]. Depending on its hydration status, it shows a hydrophilic or hydrophobic behavior: mucilage is hydrophilic in the water-saturated, swollen state, whereas it becomes hydrophobic after drying and forms hydrophobic structures on the surface of the soil particles [15,23,24]. The morphology of these structures can vary, depending on their type and concentration, from isolated local spots to networks extending over larger areas of the soil surface [12] and into the pore space [15,25,26].
For the investigation of the hydraulic properties of porous systems, such as the rhizosphere, computational fluid dynamics (CFD) is becoming an economical and effective complement or even alternative to the use of often very expensive imaging techniques such as microscale computed tomography and magnetic resonance imaging. For instances, the influence of various environmental factors on water retention curves [27,28], unsaturated hydraulic conductivity [29] or microbial diversity and activity in soils [30,31] can be analyzed by selectively varying appropriate parameters. Within the framework of the application of CFD on the pore scale, frequently used modeling methods for the investigation of fluid flows in porous media are, e.g., the finite element method [32,33] and the Lattice Boltzmann Method (LBM) [34,35,36,37,38].
To model the dynamics of spontaneous infiltration into complex porous structures with heterogeneous CA distributions, LBM is used in this work for its versatility in modelling problems with multi-phase interfaces and complex pore geometries. At the pore scale, LBM is therefore used to examine water infiltration and evaporation [39,40,41] or to obtain soil water retention characteristics [27]. Among the four most common LBMs, which are the color-liquid model [42], the Shan-Chen model (SC model) [43,44], the free-energy model [45,46] and the phase field model [47]. The SC model is utilized in this study and to simulate multiphase flows in porous media due to its higher computational efficiency and simplicity compared to the other types of models [48,49]. In addition, the model is suitable for the investigation of the infiltration process in case of heterogeneous CA distributions since the surface tension can be determined and different CA can be implemented [47]. Such infiltration processes can already be successfully simulated for simple geometric structures [50] as well as for complex porous systems [49,51,52].
Despite the important role of spatial wettability distribution for the hydraulic properties of porous media, which is intensively investigated e.g., within petroleum science [51,53,54], most simulations in soil science assume homogeneous CA for simplicity [55,56]. Only a small number of studies carry out investigations with heterogeneous wettability. They show experimentally that even hydrophilization of a few particles strongly reduces capillarity [32,57]. Of these studies, only a few deal with the effects within the rhizosphere [39,58], which is why quantitative experimental data and simulations of the influence of heterogeneous wettability on water infiltration in the rhizosphere are correspondingly rare. The aim of this study is to close this gap using a geometrically simplified model of a sandy soil and to demonstrate its practical relevance via pore scale simulations of a real sandy soil under mucilage influence.
The remainder of the paper is structured as follows. In Section 2 first the classical theory for effective contact angles for porous media with heterogenous distribution of wettability is recalled. This is followed by the details on the LBM model used in this study to investigate the infiltration at pore-scale and model validation against analytical solution. Finally, two- and three-dimensional simulation cases used in this study to systematically evaluate effect of heterogeneous wettability distributions are described. The results and discussion from the simulations are presented in Section 3. Section 3.1 deals with the dependence of infiltration on heterogeneous wettability distributions, which has been systematically investigated in over 400 simulations using a highly simplified sandy soil system. These results are compared and evaluated with those from corresponding simulations with effective contact angles. Section 3.2 describes the influence of the different wettability distributions on the water repellency of the simplified sandy soil investigated. Finally, in Section 3.3 the influence of mucilage-induced heterogeneous wettability distributions within the rhizosphere is demonstrated using a sandy soil. The conclusions from this study are finally presented in Section 4.

2. Method

This section details the classical theory for effective contact angle for heterogenous wettability distribution in the porous media and methodology for numerical simulation. In Section 2.1 the theory of the effective contact angle is first presented, as it is often used to obtain an effective homogeneous CA from heterogeneous wettability distributions. Section 2.2 describes the LBM model used to simulate infiltration process at pore-scale considering heterogenous distribution of contact angle. The model is validated against analytical solution and results are presented in Section 2.3. Subsequently, Section 2.4 describes the simulation setups used to investigate the influence of the spatial heterogeneous wettability on water infiltration at pore scale. As a first step a highly simplified, two-dimensional pore system of a sandy soil is used. Distribution of heterogenous wettability is changed systematically while keeping porous media same and in all 400 simulation variations are considered. For selected variations of heterogeneous wettability three dimensional simulations are carried out again on simplified porous media geometry to corroborate results simulated in two dimensions. Finally, Section 2.5 describes the illustrative two-dimensional simulation setup used for simulations on dry, mucilage amended sandy soil with pore-structure obtained using scanning electron microscope (SEM) to demonstrate the effect of heterogeneous wetting distributions in a rhizosphere system.

2.1. Theory of Effective Contact Angles

Often the real heterogeneous wettability distribution in porous media is represented via an effective contact angle (CA) obtained via analytical theories where the real spatial distribution is usually neglected, and averaged values are considered. One way to obtain such effective contact angles φ e f f is to use the Cassie equation [59] which was developed for flat surfaces:
cos φ e f f = A 1 cos φ 1 + A 2 cos φ 2 ,
where A 1 and A 2 are the fractions of the surface with a CA of φ 1 or φ 2 , respectively.
However, it should be noted that these theories do not take into account the spatial arrangement of the different wettabilities, but only macroscopically the area fractions of the existing wettabilities. So, to evaluate whether such effective CA can be used for simulations of infiltration processes in porous systems with heterogeneous CA distribution, φ e f f was calculated accordingly and simulations with the corresponding homogeneous CA φ e f f were carried out.

2.2. Multiphase Lattice Boltzmann Model

In this study standard pseudopotential-based LBM model for single component two phase flow is used to simulate the infiltration into two-dimensional porous media with heterogeneous CA distributions. The model has been implemented using an open-source lattice Boltzmann method framework Yantra [60,61]. In LBM model fluid is described through the evolution of particle distribution function as follows
f i x + e i t , t + t = f i x , t + t τ N S f i e q x , t f i x , t ,     i 0 , 9  
where f i represents the particle distribution function along direction i in the velocity space,   t is the time step, τ N S is relaxation time and is set to 1 for all simulation, and the set of discrete velocities e i depends on the type of lattice chosen. In this study it is the D2Q9 and D3Q19 lattice [47]. f i e q is the equilibrium distribution function given by
f i e q x , t = w i ρ 1 + e i u e s 2 + e i u 2 e s 4 u 2 e s 2 .  
For the D2Q9 lattice i 0 , 8 and the weighting factors w i = 1 / 9 for i = 1 4 , w i = 1 / 36 for i = 5 8 and w 0 = 4 / 9 and e s = e / 3 . For the D3Q19 lattice i 0 , 18 and w i = 1 / 18 for i = 1 4 , 9 , 14 , w i = 1 / 36 for i = 5 8 , 10 13 , 15 18 and w 0 = 1 / 3 and e s = e / 3 . The macroscopic quantities such as density ρ and velocity u from particle’s distribution functions as follows
ρ = i f i ,         u = 1 ρ i f i e i σ ,  
and the macroscopic kinematic viscosity ( ν ) is related to the relaxation time as
ν = e s 2 τ N S   t 2 .  
The body force F b is implemented by adding τ F b ρ to the velocity computed from Equation (4) before computing the equilibrium distribution function. For the present model the body forces are given as
F b = F i n t + F g .  
F i n t is the interparticle force accounting for non-ideal pressure during two component flow which leads to phase change and F g accounts for gravity. F i n t is given as
F i n t x , t = G   ψ x , t i w i ψ x + e i Δ t , t e i ,  
where G = 120 is the interaction strength. ψ is the interaction potential given as per the Shan-Chen model [47]:
ψ ρ = ψ 0 e ρ 0 ρ .  
In above equation ψ 0 and ρ 0 are arbitrary constant which controls the shape of equation of state. For this study, ψ 0 = 4   mu 1 2   ts   lu 1 2 and ρ 0 = 200   mu   lu 3 are chosen in such a way that the phase separation can take place with G = 120 [62,63], where mu is mass unit, lu length unit and ts time steps. For these set of parameters, the density of vapor ρ v = 85.86   mu   lu 3 and that of water ρ W = 524.98   mu   lu 3 . The density ratio between vapor and water phases is 6.11, which is significantly lower than the expected density ratio between water and air. However, the low-density ratio helps to mitigate spurious current. Moreover in this study, the mismatch between density ratio is not an issue, as flow in porous media investigated here is under a low capillary number and effects due to inertia and gravity are negligible [62]. This is further confirmed by the agreement of benchmark simulations with the Washburn equation in Section 2.3. The resulting non-ideal pressure can be computed as
P e o s = c s 2 ρ + c s 2 G 2 ψ 2 ρ .
All simulations reported in this study are calculated in lattice units and then converted into physical units by a transformation adapted to the physical system (Appendix A).
Different CA in the presence of solid surfaces can be achieved by adjusting the pseudo density of solid phase in the simulation at the positions of the wall in the range between vapor and liquid density. For solid state density values close to the water (fluid) density, a hydrophilic surface can be achieved, whereas values close to the vapor density leads to a hydrophobic surface (Figure 1).
The dependence of the equilibrium contact angle φ on the pseudo wall density ρ w a l l (Figure 2) was determined by modelling the droplet spreading on a flat surface with homogeneous wall density. Spherical cape method [64] was employed for measuring CA from equilibrium shape.

2.3. Model Validation

To validate pseudopotential LBM model implementation we test its ability to simulate a moving contact line problem against the capillary intrusion test [65]. Here, the course of the contact line between parallel plates is simulated for a two-dimensional system. Infiltration is controlled by capillary forces and the viscosity of the penetrating liquid. Under consideration of the capillary geometry and neglecting inertia, gravity, and viscosity of the gas, this leads to the established Washburn equation. It describes the following relationship for the contact line progression within the capillary [66,67]:
d l d t = γ   d cos φ   6 μ l ,  
where l is the position of the contact line, d is the plate spacing and μ is the dynamic viscosity of the liquid. For a homogeneous coating of the panels with one CA , the analytical solution of Equation (10) is:
l = γ   d cos φ   3 μ t .
A system of 7 × 0.2   mm 2 with a conversion factor for length C L = 10 μ m lu is used for simulation. The parallel plates with a length of 5   mm are arranged in a way that they are at a distance of 0.1   mm from each other and 1   mm each in the x-direction from the side walls (Figure 3). Periodic boundary conditions are used in the y-direction where there are no parallel plates. The parallel plates are implemented as solid obstacles with bounce-back boundary conditions. At the left side of the system, a Dirichlet condition is set to the density of water, representing an infinite water reservoir. Whereas at the right edge the Neumann condition ( u = 0 ) prevents a drop formation at the right end of the plates.
As initial condition of the simulation the liquid is placed on the left side of the plates. After a short system settling of 7   μ s , while the liquid meniscus slightly penetrated between the plates up to point x 0 close to the entry of the capillary, the speed u was set back to zero in the entire system. Subsequently, the progression of the meniscus in relation to the reference point x 0 was measured as a function of time.
The results for homogeneous CA of 0 ° and 60 ° are shown in Figure 4 and indicate that the penetration curves are basically divided into two different stages when comparing Washburn and LBM results. In the first stage up to about 0.08   ms , the penetration velocity of the LBM is significantly lower and therefore more reflects the actual course of the initial phase [68]. In the second stage, the LBM results approximately describe the Washburn process but the deviation increases over time.
To validate this LBM model also for spatial heterogeneous CA , the same system as before was simulated with alternating coating of the parallel surfaces, i.e., two different CA follow each other at an even spatial distance in flow direction. (Figure 5). The CA pairs used are 0° and 30° as well as 0° and 60° with an equidistant strip width of 0.2   mm .
With regard to the results of the penetration at alternating CA , two phases can be identified, as in the case of homogenous coating. In the first phase ( t < 0.085   ms ), the average speed derived from the LBM simulations is lower while in the second phase it is higher than the numerical solution of the Washburn equation. The spatial distribution of the alternating CA is clearly reflected in the different penetration velocities and matches them, even though the temporal fit may vary in sections due to the time-dependent velocity difference mentioned above.

2.4. Infiltration into Porous Media Depending on Contact Angle Distribution

To examine the influence of CA distributions on infiltration into porous media, at first, the infiltration into geometrically highly simplified porous systems is simulated using the previously described and validated program. The 1.17 × 0.17   mm 2 systems employed here with a conversion factor for length C L = 1.7 μ m lu are identical to those utilized for the capillary intrusion test in terms of boundary conditions and periodicity. The main difference lies in the fact that the parallel plates are replaced by a porous system consisting of circular particles positioned in a regular way (Figure 6). The porosity is adjusted to resemble that of sand with Φ = 43.4   % and the unit conversion was based on a sand particle diameter of d = 0.1   mm .
The CA coating of the particles is applied according to three different coating patterns (Table 1, Figure 7):
CA of 0°, 20°, 30°, 40°, 50°, 55°, 60°, 70°, 80 and 90° are used to investigate the dynamics of the infiltration process. Two different wettability patterns are investigated in the cases of the body-throat coating pattern: in the body-phobic-case, the pore throat exhibits a hydrophilic CA of 0° and the pore body exhibits various reduced wettabilities (Figure 7b). In the throat-phobic-case, the wettabilities of the pore throat and body are reversed. In addition to the CA , the influence of the size of the surface with reduced wettability is also investigated for all body-throat simulations and is varied in the form of a surface angle α from 5° up to 60° in 5° steps. For the body-phobic coating, α describes the strip width with reduced wettability located in the pore body, while for the throat-phobic it is related to the strip width with reduced wettability in the pore throat. For coating pattern stripes, hydrophilic stripes are always followed by a stripe with reduced wettability whereby the stripe width is defined by the angle α of 5°, 10° or 15° (Figure 7c).
To compare the influence of the coating pattern on dynamic infiltration, the average front velocity v ¯ of the fluid was determined by measuring the time until it penetrated 0.58   mm into the porous medium. For each coating pattern, the critical contact angle (CCA) at which no infiltration takes place was determined by additional simulations to an accuracy of 1 ° .
In addition, the infiltration is simulated simplified, three-dimensional porous medium 1.17 × 0.17 × 0.17   mm 3 (Figure 8). Boundary conditions and periodicity are now adjusted to three dimensions and they are identical for the y- and z-direction. The CA distribution from the two-dimensional throat-phobic-case is now extended to three dimensions. For this purpose, the wettability is reduced at the surface stripes of the spherical particles, which have the smallest distance between neighboring spheres and are oriented vertically to the direction of infiltration (x-direction). Also experimentally it has been observed that mucilage during drying retreats into the locations of smallest distances between particles [15,25]. For comparability with infiltration into soil from Section 2.5, a CA of 100° was chosen for these surfaces, while the rest has a CA of 20°. Infiltration scenarios were carried out for surface strip widths α of 5°, 10° or 15°.

2.5. Infiltration into Soil

To demonstrate the dependence of water infiltration on the heterogeneous nature of CA distribution of soils, the infiltration process was simulated based on the example of a dry, sandy soil (0.125–0.2 mm) amended with maize (Zea mays L.) mucilage at a content of 8   mg   g 1 ( mg dry mucilage per g of dry soil). To resolve the spatial distribution of mucilage induced by drying in soil, the sand was amended with hydrated mucilage and dried by evaporation (see [25] for a detailed description). The distribution of dry mucilage structures created in the process of soil drying was resolved using synchrotron-based X-ray tomographic microscopy (SRXTM). Mucilage structures were visible as two-dimensional surfaces connecting multiple pores across the soil domain [25]. The origin of two-dimensional structures on soil particles is indicated in green on the example displayed below (Figure 9; courtesy of M. Zarebanadkouki, University of Bayreuth, and P. Benard, A. Carminati, ETH Zurich).
The simulations were carried out for this two-dimensional partial section of the sandy soil (Figure 9) with a lattice of 0.6 × 0.48   mm 2 and a conversion factor for length C L = 0.33 μ m lu . The boundary conditions as well as the periodicity of the simulations correspond to those of the simplified porous medium (Figure 6), whereby two thin parallel plates were additionally added above and below the soil section. The plates have a CA of 90° to ensure that, at the top and bottom boundaries, the curvature of the water front behaves as if the system was mirrored in the y-direction. The CA of the areas covered with mucilage is set to 100 ° , which corresponds to observed values at high maize mucilage surface concentrations [23]. These values can be expected locally as the retreat of liquid during drying leads to a high concentration of mucilage in this regions [25,26]. The CA of the remaining area of soil particle surfaces was set to 20 ° , a typical value for untreated soil mineral surfaces. The area covered with mucilage is 5.2% of the total surface area of the particles. For this distribution, the Cassie equation predicts an effective CA of 28°.
Two simulations were performed: one with the heterogeneous distributions of CA as described above and one assuming a homogeneous effective CA of 28° according to the Cassie equation. To exclude a filling of the soil section due to artificial condensation (DeMaio2011), the determination of the equilibrium is here based exclusively on the fluid, which is in direct contact with the water reservoir.

3. Results and Discussion

3.1. Dynamic: Simplified Porous System

3.1.1. Homogeneous Coating Pattern

For the infiltration process of the homogeneous coating pattern case, the time-dependent penetration length directly reflects the geometry of the porous system (Figure 10 and Figure 11): the infiltration front velocity v ¯ in the pore throat is significantly greater than in the pore body. This can be explained by the interplay of continuity equation which states a higher flow velocity in narrow parts of the system and the stronger capillary forces in the pore throat compared to the pore body. Note that for infiltration in a capillary of constant radius, instead, the Washburn equation (Equation (11)) predicts a lower v ¯ at smaller diameters.
The mean velocity v ¯ decreases according to the Washburn equation and the CCA is 54° (Figure 10 and Figure 12). For values above the CCA water cannot infiltrate into the system and the average velocity is 0.
Transferring this to soils in general, this means that water repellency can occur not only for hydrophobic CA   ( > 90 ° ) but also at reduced wettability of the soil particles. The reason for this is that due to the system geometry the macro-scale, effective CA depends not only on the curvature of the water surface (as with a flat surface), but also on the curvature of the particle surface [69].

3.1.2. Body-Throat Coating Patterns

Also, for the body-throat coating patterns (Figure 13), in general, v ¯ decreases with increasing CA . The small increase in throat-phobic coatings for α 20 ° (Figure 13b) are caused by discretization effects in the pore throat. For same CA , v ¯ decreases with increasing α . It should be emphasized that for body-phobic, for α 20 ° , the value of α has nearly no influence on v ¯ and also corresponds approximately to the homogeneous coating pattern (Figure 13a). On the other hand, for α < 20 °   v ¯ is directly dependent on α . For the throat-phobic pattern, in contrast, v ¯ depends strongly on α , for α 20 ° , but is more or less independent of CA for α < 20 ° (Figure 13b).
A comparison of v ¯ of body-phobic and throat-phobic (Figure 14) shows that at same α , v ¯ in the body-phobic-case is always smaller than in the throat-phobic-case. If we consider v ¯ as a measure of the infiltration rate σ of the porous system and take into account that the total infiltration rate is limited by the globally lowest infiltration rate in this system, we can see that for a composite system of body-phobic and throat-phobic with α b o d y - p h o b α t h r o a t - p h o b the effective total infiltration rate σ e f f corresponds to the infiltration rate of a system with only body-phobic coating. For α b o d y - p h o b 20 ° , σ e f f is even determined independently of α b o d y - t h r o a t of the infiltration rate of body-phobic and can be approximately calculated by a homogeneous coating with a CA equal to the CA of body-phobic.
Applied to porous media, these results show that the infiltration dynamics strongly depend on the CA distribution. Furthermore, the wettability of the pore body is here more decisive for the infiltration than the wettability of the pore throat. However, it should be noted that the two-dimensional findings cannot be directly transferred to the three-dimensional real soil system case: in three dimensions the pore body has a much larger surface compared to the one of the pore throats. Therefore, much more hydrophobic material would be required to make a 3D pore body water repellent.

3.1.3. Stripes Coating Pattern

For the stripes coatings v ¯ also decreases with increasing CA , whereby for CA 40 ° it is rather independent of α (Figure 15). For α = 10 ° and 15 ° the curve is very similar up to CA 60 ° , which could possibly be caused by the fact that for α = 10 ° the opposite stripes of two neighboring particles have the same CA , whereas for α = 15 ° these are distinct (Figure 7c).
v ¯ of the strip coating is smaller than that of the body-phobic for the same α which can be explained by the smaller fraction of the total surface with reduced wettability of the strips near the pore body (Figure 16). For α = 15 ° this is reversed, which is probably due to the different CA of the stripes facing each other.
From results of the stripes simulations, it can be concluded for porous systems in general that even at same CA area ratio a smaller characteristic pattern size of the CA distribution could increase the marco-scale effective infiltration rate.

3.1.4. Comparison with Cassie Equation

In the following, the dependence of v ¯ on CA is compared for coatings with the same effective CA but different spatial CA distribution to study the basic scope and limitations of effective CA . For this purpose, all coating samples with 1:1 ratio of hydrophilic surfaces to surfaces with reduced wettability are compared with the results of a corresponding homogeneous effective CA coating according to the Cassie equation. All results based on the homogeneous, effective CA coating are summarized in the following as Cassie curves. The comparison of the simulation results of the Cassie curve with the other coating patterns with respect to v ¯ shows that the deviations from the Cassie curve increase with increasing CA (Figure 17 and Figure 18). The Cassie curve corresponds approximately to the curves of the strip coating for CA 40 ° , whereby strong deviations are recognizable for CA > 60 ° . Additionally, the stripes coating does not converge against the Cassie equation for decreasing α , since for the stripes coating at α = 5 ° the values of v ¯ are always larger and for α = 10 ° and 15 ° always smaller than for the Cassie curve. Comparing the Cassie curve with those of body-phobic and throat-phobic, the Cassie curve is between them, whereby the difference for throat-phobic is always greater and also increases with increasing CA.
For porous media, it can be deduced that the specific spatial CA distribution pattern has a decisive influence on the infiltration process and that it is not sufficient to only take into account the relative area fraction of different CA s. Particularly in the case of relatively large CA differences as well as large α in the pore body, caution is required for the application of the Cassie equation. Moreover, the non-existence of convergence for small α against the Cassie curve generally raises the question about the reasonable application of the Cassie equation in more complex porous media, such as in the context of soils.

3.2. Critical Contact Angle (CCA)

The CCA is greater for the throat-phobic coating than for the body-phobic coating for all α (Figure 19). In the case of the body-phobic coating it changes from strongly hydrophobic 142 ° to reduced wettability 54 ° in the range of 5 ° α 25 ° . For larger α the CCA is constant 54 ° and thus corresponds to the CCA of the homogeneous coating, which is why it is reasonable to assume that the CCA is also constant 54 ° for 60 ° < α < 90 ° . Even for the throat-phobic coating the CCA initially drops sharply in the hydrophobic range (from 159 ° to 102 ° ) for 5 ° α 10 ° . For larger α it decreases more weakly and approximately constant, changing from hydrophobic to reduced wettability in the range 20 ° α 25 ° . It is assumed to decrease further for α 60 ° until it reaches CCA = 54 ° at α = 90 ° which corresponds to the homogeneous coating. The comparison of the coatings shows that for α = 5 ° the CCA of the stripes coating with 90° is much smaller than for the body-phobic and throat-phobic coating, which may be explained by the additional stripes of the stripes coating. The CCA for α 10 ° is approximately equal to that of the body-phobic coating, since the stripes in the pore body are more relevant (Figure 16).
Regarding the behavior for small α , Figure 19 supports two hypotheses: we expect that for the body-throat coating at sufficiently small α , a CCA no longer exists, since the momentum of the infiltration front allows the fluid to flow over a sufficiently small hydrophobic layer even at a CA of 180 °. For the stripes coating on the other hand, it would be conceivable that for small α the CCA converges towards a fixed value.
Transferred to porous media, the water repellency can therefore strongly depend on the spatial CA distribution. In our two-dimensional cases, the wettability of pore bodies is more decisive than the one of pore throats. In order to transfer the data to three-dimensional systems, two main cases have to be considered, which can be distinguished on the basis of the CA . If only reduced wettability is present, water repellency is largely controlled by the coating of the pore body. In the second case, if additional hydrophobic areas occur in the pore throat, these can lead to water repellency. These areas in the pore throat are especially important because they control water content dynamics under unsaturated conditions. The three-dimensional simulations (Figure 20) showing water repellency from a surface stripe width of α = 15 ° underline that even small hydrophobic surface fractions (here 17 % of the total surface) in the pore throat are sufficient to cause water repellency.

3.3. Cassie Equation in Soil

Water infiltration in sandy soil is highly dependent on the spatial distribution of wettability and not only on the specific area ratio. Simulated infiltration into a soil system covered at certain spots with mucilage (with mucilage: CA = 100 ° without: CA = 20 ° ) shows that the system is water repellent when covered partially with mucilage (Figure 21), while it is conductive with a corresponding homogeneous coating with the equivalent CA = 28 ° according to the Cassie equation. Thus, this underlines the hypothesis that for complex porous media theories like the Cassie equation that are based on an effective CA and do not consider the specific distribution may be unsuitable for estimating infiltration processes and that even small areas covered with dry mucilage (here 5% of the total surface) can have a considerable influence on water repellency.

4. Conclusions

In our study we investigated the influence of spatially heterogeneous wettability distributions on infiltration processes in simplified and real two-dimensional sandy soil systems at pore-scale. We not only confirmed that water repellency in porous media such as soil can occur for non-hydrophobic CA [70] but also we could show that water repellency as well as water infiltration dynamics are strongly influenced by the spatial distribution of wettability. We showed that simulations based on the effective contact angle theory such as the Cassie equation can predict dynamics that are far from the results obtained when considering location and patterns of local contact angle distributions.
At certain conditions, such as big CA variations, the Cassie equation is insufficient for calculating effective contact angles for heterogeneous CA distributions in porous media [71,72]. For two-dimensional systems we could show that, the infiltration and water repellency is more sensitive to increases in contact angle and area of coated surfaces in pore bodies than to increases in pore throats. Hydrophobic locations in the pore throat, but also in the pore body, are each sufficient to cause water repellency. Indeed, we could show for a two-dimensional real soil system that already covering certain pore throats which corresponds to 5% of the total surface with a hydrophobic layer is sufficient to induce water repellency. In the three-dimensional case, pore throats comprise far less soil surface compared to pore bodies. This means that there even smaller fractions of the soil surface in pore throats covered with a hydrophobic material can render soil water repellency, especially in unsaturated cases. For instance, in the rhizosphere during drying, mucilage retreats to the narrower parts of the pore throats and finally dries there on the soil surface [15], covering the pore throats with a hydrophobic layer of a contact angle of sometimes around 100°. We have shown, that even very tiny amounts of such a layer in the proper location can induce water repellency. Since pore throats are the last parts of a soil (before surfaces) to dry and the dependence of CA on water content [73], it is also conceivable that a critical water content exists below which the respective soil, and in particular the rhizosphere, is water repellent. With regard to the application to natural soil, it should be noted that the described effects can vary considerably depending on particle geometry, roughness and contact angle.
Our study shows that averaged contact angle parametrizations are not appropriate and a proper knowledge about local contact angle patterns needs to be considered to advance our understanding of the effect of microscale properties on large scale hydraulic dynamics and to include the microscale effect of root exudation into larger root water uptake models [74].
Finally, it should be noted that this study focuses on investigating the influence of the wettability distribution when the pore structure properties of the porous medium are constant. The conclusions derived from this study can apply to all porous media in general. However, further simulations are needed to quantitatively evaluate the additional influence of the heterogeneity of the porous media in combination with the wettability distribution.

Author Contributions

J.B.: extended the simulation software, created and ran simulations, analyzed data, developed conceptual model, wrote original draft preparation. E.K.: supervision, conceptualization, writing—reviewing and editing. R.A.P.: extended the simulation software, conceptualization, writing—reviewing and editing. P.B.: provided SRXTM images, writing—reviewing and editing. A.L.: determined mucilage distribution on the basis of the SRXTM images. A.H.: conceptualization. All authors have read and agreed to the published version of the manuscript.

Funding

J.B. was funded by the MUCI-WETT-PATT Project in the frame of the priority program PP2089 “Rhizosphere Spatiotemporal Organization” (403668613). The position of P.B. was funded by the German Research Foundation DFG under the project number 403640522 within the project “Liquid Rhizosphere” of the priority program 2089 “Rhizosphere spatiotemporal organization—a key to rhizosphere functions”. The positions of A.H. and E.K. have been funded by the German Research Foundation under Germany’s Excellence Strategy, EXC-2070—390732324—PhenoRob. A.L. was funded within the framework of the priority program 2089, funded by the DFG, project 403660839.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The general data are included in the article. Additional data are available on request.

Acknowledgments

We acknowledge the German Research Foundation priority program PP2089 for funding this project.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

CAcontact angle
CCAcritical contact angle
CFDcomputational fluid dynamics
LBMLattice Boltzmann Method
lu length unit
mu mass unit
SCShan-Chen
SEM scanning electron microscope
SOMsoil organic matter
SRXTMsynchrotron-based X-ray tomographic microscopy
SWR soil water repellency
ts time steps
WB:Washburn
Symbols:
A surface fraction
c grid spacing
C L conversion factor of length
C M conversion factors of mass
C T conversion factors of time
C γ conversion factors of surface tension
C μ conversion factors of dynamic viscosity
d plate spacing or particle diameter
e i discrete velocities along direction i
F b body force
F g gravitational force
F i n t interparticle force
f i   particle distribution function along direction i
f i e q equilibrium particle distribution function along direction i
G interaction strength
l position of the contact line
L physical length
L 2 error norm
N number of lattice nodes
P pressure
R 2 correlation coefficient
t time
u macroscopic velocity
v ¯ infiltration front velocity
w i weighting factor along direction i
x penetration length
x f i n e s t penetration length of finest grid
x 0 penetrated length after short system settling
α stripe width
α b o d y - p h o b stripe width of the body of a body-phobic coating
α t h r o a t - p h o b stripe width of the body of a body-phobic coating
γ surface tension
t time step
Φ porosity
μ dynamic viscosity
ν kinematic viscosity
ρ macroscopic density
ρ 0 constant of interaction potential function
ρ w a l l pseudo density of the wall
σ e f f effective total infiltration rate
τ NS relaxation time
φ contact angle
φ e f f effective contact angle
ψ interaction potential
ψ 0 constant of interaction potential function

Appendix A. Unit Conversion

The system variables represented in the LBM in grid units (lb) can be converted into physical units (phy) with the help of conversion factors. Starting from the basic conversion factors for length C L , time C T and mass C M , all quantities derived from these can be directly calculated. The conversion for infiltration processes can be done using dimensionless numbers such as the Bond or Reynolds number which can be adjusted for our simulations with help of system length, dynamic viscosity and surface tension. The length conversion factor is defined as the ratio between the physical length L of the computational domain and the corresponding number of lattice nodes N in the grid:
C L = L p h y N .
With the help of the conversion factors of the dynamic viscosity C μ of fluid
C μ = C ρ   C ν = C M C L 3 C L 2 T = C M C L   C T ,
and the surface tension C γ
C γ = C M C T 2 ,
the conversion factors for time and mass can be calculated by reshaping and inserting Equations (A2) and (A3):
C T = C μ C γ C L ,  
C M = C μ 2 C γ C L 2 .
The kinetic viscosity in lattice units ν l b is determined by Equation (5) with τ N S = 1 , which leads to ν l b = 1 6   lu 2   ts 1 for water and vapor. Using the densities of vapor ρ v = 85.86   mu   lu 3 and fluid ρ W = 524.98   mu   lu 3 , the dynamic viscosity of vapor is thus μ l b , v = 14.31 mu lu   ts and that of the fluid   μ l b , w = 87.5 mu lu   ts .
The calculation of the surface tension γ in this multi-phase LB-model in lattice units is done based on the two-dimensional Young-Laplace equation
P = γ r ,  
by determining the pressure difference P in bubbles or droplets of different radii r . Since the interface transition between vapor and liquid is continuous, the arithmetic mean of the equilibrium vapor and liquid density was chosen as the limiting density for phase differentiation ρ I = 305.42   mu   lu 3 . A linear fitting of the points on the curve P ~ 1 / r results in a surface tension of γ = 13.85   mu   ts 2 with correlation coefficient R 2 = 0.9988 (Figure A1).
Figure A1. Correlation between pressure difference ( p ) between inside and outside the bubble and reciprocal of its equilibrium radius r 1 as well as the resulting surface tension γ .
Figure A1. Correlation between pressure difference ( p ) between inside and outside the bubble and reciprocal of its equilibrium radius r 1 as well as the resulting surface tension γ .
Water 14 01110 g0a1

Appendix B. Grid Convergence

To check the grid convergence, the same simulation setup as for the capillary intrusion test with homogeneous hydrophilic wetting was used. The grid spacing is defined by c 1 = 1 × 10 5   mm , c 2 = 5 × 10 6   mm , c 3 = 7 × 10 7   mm and c 4 = 1 × 10 7   mm .
The relative error (Figure A2) between the grid with the finest grid spacing c 4 = 1 × 10 7   mm and the other coarser grids is represented in form of the L 2 error norm and can be calculated as follows
L 2 = x x f i n e s t 2 x f i n e s t 2
with x and x f i n e s t as the time dependent penetration length of the coarser grid and finest grid respectively (Figure A2b).
Figure A2. Illustration of convergence with respect to grid spacing using the simulation setup for the capillary intrusion test of water into parallel plates. The figure presents the effect of the grid spacing on (a) the curves of time-dependent penetration length and (b) the corresponding L 2 diagram with a convergence order of 1.
Figure A2. Illustration of convergence with respect to grid spacing using the simulation setup for the capillary intrusion test of water into parallel plates. The figure presents the effect of the grid spacing on (a) the curves of time-dependent penetration length and (b) the corresponding L 2 diagram with a convergence order of 1.
Water 14 01110 g0a2

References

  1. Vereecken, H.; Schnepf, A.; Hopmans, J.W.; Javaux, M.; Or, D.; Roose, T.; Vanderborght, J.; Young, M.H.; Amelung, W.; Aitkenhead, M.; et al. Modeling Soil Processes: Review, Key Challenges, and New Perspectives. Vadose Zone J. 2016, 15, vzj2015.09.0131. [Google Scholar] [CrossRef] [Green Version]
  2. Vogel, H.-J.; Roth, K. Moving through Scales of Flow and Transport in Soil. J. Hydrol. 2003, 272, 95–106. [Google Scholar] [CrossRef]
  3. Doerr, S.H.; Shakesby, R.A.; Walsh, R.P.D. Soil Water Repellency: Its Causes, Characteristics and Hydro-Geomorphological Significance. Earth-Sci. Rev. 2000, 51, 33–65. [Google Scholar] [CrossRef]
  4. Young, T. III. An Essay on the Cohesion of Fluids. Philos. Trans. R. Soc. Lond. 1805, 95, 65–87. [Google Scholar] [CrossRef]
  5. Tschapek, M. Criteria for Determining the hydrophilicity-hydrophobicity of Soils. Z. Pflanz. Bodenkd. 1984, 147, 137–149. [Google Scholar] [CrossRef]
  6. Ma’Shum, M.; Tate, M.E.; Jones, G.P.; Oades, J.M. Extraction and Characterization of Water-Repellent Materials from Australian Soils. J. Soil Sci. 1988, 39, 99–110. [Google Scholar] [CrossRef]
  7. Woche, S.K.; Goebel, M.-O.; Mikutta, R.; Schurig, C.; Kaestner, M.; Guggenberger, G.; Bachmann, J. Soil Wettability Can Be Explained by the Chemical Composition of Particle Interfaces—An XPS Study. Sci. Rep. 2017, 7, 42877. [Google Scholar] [CrossRef]
  8. Cerdà, A.; Doerr, S.H. Soil Wettability, Runoff and Erodibility of Major Dry-Mediterranean Land Use Types on Calcareous Soils. Hydrol. Process. 2007, 21, 2325–2336. [Google Scholar] [CrossRef]
  9. Goebel, M.-O.; Bachmann, J.; Woche, S.K. Modified Technique to Assess the Wettability of Soil Aggregates: Comparison with Contact Angles Measured on Crushed Aggregates and Bulk Soil. Eur. J. Soil Sci. 2008, 59, 1241–1252. [Google Scholar] [CrossRef]
  10. Madsen, M.D.; Zvirzdin, D.L.; Petersen, S.L.; Hopkins, B.G.; Roundy, B.A.; Chandler, D.G. Soil Water Repellency within a Burned Piñon-Juniper Woodland: Spatial Distribution, Severity, and Ecohydrologic Implications. Soil Sci. Soc. Am. J. 2011, 75, 1543–1553. [Google Scholar] [CrossRef] [Green Version]
  11. Bachmann, J.; Goebel, M.-O.; Woche, S.K. Small-Scale Contact Angle Mapping on Undisturbed Soil Surfaces. J. Hydrol. Hydromech. 2013, 61, 3–8. [Google Scholar] [CrossRef] [Green Version]
  12. Kaltenbach, R.; Diehl, D.; Schaumann, G.E. Links between Nanoscale and Macroscale Surface Properties of Natural Root Mucilage Studied by Atomic Force Microscopy and Contact Angle. J. Colloid Interface Sci. 2018, 516, 446–455. [Google Scholar] [CrossRef]
  13. Gregory, P.J. Roots, Rhizosphere and Soil: The Route to a Better Understanding of Soil Science? Eur. J. Soil Sci. 2006, 57, 2–12. [Google Scholar] [CrossRef]
  14. Bais, H.P.; Weir, T.L.; Perry, L.G.; Gilroy, S.; Vivanco, J.M. The role of root exudates in rhizosphere interactions with plants and other organisms. Annu. Rev. Plant Biol. 2006, 57, 233–266. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  15. Carminati, A.; Benard, P.; Ahmed, M.A.; Zarebanadkouki, M. Liquid Bridges at the Root-Soil Interface. Plant Soil 2017, 417, 1–15. [Google Scholar] [CrossRef]
  16. Hinsinger, P.; Bengough, A.G.; Vetterlein, D.; Young, I.M. Rhizosphere: Biophysics, Biogeochemistry and Ecological Relevance. Plant Soil 2009, 321, 117–152. [Google Scholar] [CrossRef]
  17. Whalley, W.R.; Riseley, B.; Leeds-Harrison, P.B.; Bird, N.R.A.; Leech, P.K.; Adderley, W.P. Structural Differences between Bulk and Rhizosphere Soil. Eur. J. Soil Sci. 2005, 56, 353–360. [Google Scholar] [CrossRef]
  18. Carminati, A.; Moradi, A.B.; Vetterlein, D.; Vontobel, P.; Lehmann, E.; Weller, U.; Vogel, H.-J.; Oswald, S.E. Dynamics of Soil Water Content in the Rhizosphere. Plant Soil 2010, 332, 163–176. [Google Scholar] [CrossRef]
  19. Kroener, E.; Holz, M.; Zarebanadkouki, M.; Ahmed, M.; Carminati, A. Effects of Mucilage on Rhizosphere Hydraulic Functions Depend on Soil Particle Size. Vadose Zone J. 2018, 17, 170056. [Google Scholar] [CrossRef] [Green Version]
  20. Kroener, E.; Zarebanadkouki, M.; Kaestner, A.; Carminati, A. Nonequilibrium Water Dynamics in the Rhizosphere: How Mucilage Affects Water Flow in Soils. Water Resour. Res. 2014, 50, 6479–6495. [Google Scholar] [CrossRef] [Green Version]
  21. Muñoz, L.A.; Cobos, A.; Diaz, O.; Aguilera, J.M. Chia Seeds: Microstructure, Mucilage Extraction and Hydration. J. Food Eng. 2012, 108, 216–224. [Google Scholar] [CrossRef]
  22. Read, D.B.; Gregory, P.J.; Bell, A.E. Physical Properties of Axenic Maize Root Mucilage. Plant Soil 1999, 211, 87–91. [Google Scholar] [CrossRef]
  23. Ahmed, M.A.; Kroener, E.; Benard, P.; Zarebanadkouki, M.; Kaestner, A.; Carminati, A. Drying of Mucilage Causes Water Repellency in the Rhizosphere of Maize: Measurements and Modelling. Plant Soil 2016, 407, 161–171. [Google Scholar] [CrossRef]
  24. Moradi, A.B.; Carminati, A.; Lamparter, A.; Woche, S.K.; Bachmann, J.; Vetterlein, D.; Vogel, H.-J.; Oswald, S.E. Is the Rhizosphere Temporarily Water Repellent? Vadose Zone J. 2012, 11, vzj2011.0120. [Google Scholar] [CrossRef]
  25. Benard, P.; Zarebanadkouki, M.; Brax, M.; Kaltenbach, R.; Jerjen, I.; Marone, F.; Couradeau, E.; Felde, V.J.M.N.L.; Kaestner, A.; Carminati, A. Microhydrological Niches in Soils: How Mucilage and EPS Alter the Biophysical Properties of the Rhizosphere and Other Biological Hotspots. Vadose Zone J. 2019, 18, 180211. [Google Scholar] [CrossRef] [Green Version]
  26. Benard, P.; Zarebanadkouki, M.; Hedwig, C.; Holz, M.; Ahmed, M.A.; Carminati, A. Pore-Scale Distribution of Mucilage Affecting Water Repellency in the Rhizosphere. Vadose Zone J. 2018, 17, 170013. [Google Scholar] [CrossRef] [Green Version]
  27. Ahrenholz, B.; Tölke, J.; Lehmann, P.; Peters, A.; Kaestner, A.; Krafczyk, M.; Durner, W. Prediction of Capillary Hysteresis in a Porous Material Using Lattice-Boltzmann Methods and Comparison to Experimental Data and a Morphological Pore Network Model. Adv. Water Resour. 2008, 31, 1151–1173. [Google Scholar] [CrossRef]
  28. Garg, A.; Garg, A.; Tai, K.; Barontini, S.; Stokes, A. A Computational Intelligence-Based Genetic Programming Approach for the Simulation of Soil Water Retention Curves. Transp. Porous Media 2014, 103, 497–513. [Google Scholar] [CrossRef]
  29. Daly, K.R.; Mooney, S.J.; Bennett, M.J.; Crout, N.M.J.; Roose, T.; Tracy, S.R. Assessing the Influence of the Rhizosphere on Soil Hydraulic Properties Using X-ray Computed Tomography and Numerical Modelling. J. Exp. Bot. 2015, 66, 2305–2314. [Google Scholar] [CrossRef] [Green Version]
  30. Baveye, P.C.; Otten, W.; Kravchenko, A.; Balseiro-Romero, M.; Beckers, É.; Chalhoub, M.; Darnault, C.; Eickhorst, T.; Garnier, P.; Hapca, S.; et al. Emergent Properties of Microbial Activity in Heterogeneous Soil Microenvironments: Different Research Approaches Are Slowly Converging, Yet Major Challenges Remain. Front. Microbiol. 2018, 9, 1929. [Google Scholar] [CrossRef] [PubMed]
  31. Portell, X.; Pot, V.; Garnier, P.; Otten, W.; Baveye, P.C. Microscale Heterogeneity of the Spatial Distribution of Organic Matter Can Promote Bacterial Biodiversity in Soils: Insights from Computer Simulations. Front. Microbiol. 2018, 9, 1583. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  32. Deurer, M.; Bachmann, J. Modeling Water Movement in Heterogeneous Water-Repellent Soil: 2. A Conceptual Numerical Simulation. Vadose Zone J. 2007, 6, 446–457. [Google Scholar] [CrossRef]
  33. Jafari, I.; Masihi, M.; Nasiri Zarandi, M. Numerical Simulation of Counter-Current Spontaneous Imbibition in Water-Wet Fractured Porous Media: Influences of Water Injection Velocity, Fracture Aperture, and Grains Geometry. Phys. Fluids 2017, 29, 113305. [Google Scholar] [CrossRef]
  34. Aljasmi, A.; Sahimi, M. Speeding-up Image-Based Simulation of Two-Phase Flow in Porous Media with Lattice-Boltzmann Method Using Three-Dimensional Curvelet Transforms. Phys. Fluids 2021, 33, 113313. [Google Scholar] [CrossRef]
  35. Lin, W.; Xiong, S.; Liu, Y.; He, Y.; Chu, S.; Liu, S. Spontaneous Imbibition in Tight Porous Media with Different Wettability: Pore-Scale Simulation. Phys. Fluids 2021, 33, 032013. [Google Scholar] [CrossRef]
  36. Seta, T.; Takegoshi, E.; Okui, K. Lattice Boltzmann Simulation of Natural Convection in Porous Media. Math. Comput. Simul. 2006, 72, 195–200. [Google Scholar] [CrossRef]
  37. Guo, Z.; Zhao, T.S. Lattice Boltzmann Model for Incompressible Flows through Porous Media. Phys. Rev. E 2002, 66, 036304. [Google Scholar] [CrossRef]
  38. Abouei Mehrizi, A.; Farhadi, M.; Sedighi, K.; Latif Aghili, A. Lattice Boltzmann Simulation of Heat Transfer Enhancement in a Cold Plate Using Porous Medium. J. Heat Transf. 2013, 135, 111006. [Google Scholar] [CrossRef]
  39. Cruz, B.C.; Furrer, J.M.; Guo, Y.-S.; Dougherty, D.; Hinestroza, H.F.; Hernandez, J.S.; Gage, D.J.; Cho, Y.K.; Shor, L.M. Pore-Scale Water Dynamics during Drying and the Impacts of Structure and Surface Wettability: Pore-scale water dynamics during drying. Water Resour. Res. 2017, 53, 5585–5600. [Google Scholar] [CrossRef]
  40. Li, Y.; LeBoeuf, E.J.; Basu, P.K.; Mahadevan, S. Stochastic Modeling of the Permeability of Randomly Generated Porous Media. Adv. Water Resour. 2005, 28, 835–844. [Google Scholar] [CrossRef] [Green Version]
  41. Sukop, M.C.; Or, D. Invasion Percolation of Single Component, Multiphase Fluids with Lattice Boltzmann Models. Phys. B Condens. Matter 2003, 338, 298–303. [Google Scholar] [CrossRef]
  42. Grunau, D.; Chen, S.; Eggert, K. A Lattice Boltzmann Model for Multiphase Fluid Flows. Phys. Fluids Fluid Dyn. 1993, 5, 2557–2562. [Google Scholar] [CrossRef] [Green Version]
  43. Shan, X.; Chen, H. Simulation of Nonideal Gases and Liquid-Gas Phase Transitions by the Lattice Boltzmann Equation. Phys. Rev. E 1994, 49, 2941–2948. [Google Scholar] [CrossRef] [Green Version]
  44. Shan, X.; Chen, H. Lattice Boltzmann Model for Simulating Flows with Multiple Phases and Components. Phys. Rev. E 1993, 47, 1815–1819. [Google Scholar] [CrossRef] [Green Version]
  45. Swift, M.R.; Orlandini, E.; Osborn, W.R.; Yeomans, J.M. Lattice Boltzmann Simulations of Liquid-Gas and Binary Fluid Systems. Phys. Rev. E 1996, 54, 5041–5052. [Google Scholar] [CrossRef]
  46. Swift, M.R.; Osborn, W.R.; Yeomans, J.M. Lattice Boltzmann Simulation of Nonideal Fluids. Phys. Rev. Lett. 1995, 75, 830–833. [Google Scholar] [CrossRef] [Green Version]
  47. Huang, H.; Sukop, M.C.; Lu, X.-Y. Multiphase Lattice Boltzmann Methods: Theory and Application; John Wiley and Sons, Inc.: Chichester, UK, 2015; ISBN 978-1-118-97133-8. [Google Scholar]
  48. Li, Q.; Luo, K.H.; Kang, Q.J.; He, Y.L.; Chen, Q.; Liu, Q. Lattice Boltzmann Methods for Multiphase Flow and Phase-Change Heat Transfer. Prog. Energy Combust. Sci. 2016, 52, 62–105. [Google Scholar] [CrossRef] [Green Version]
  49. Chen, L.; Kang, Q.; Mu, Y.; He, Y.-L.; Tao, W.-Q. A Critical Review of the Pseudopotential Multiphase Lattice Boltzmann Model: Methods and Applications. Int. J. Heat Mass Transf. 2014, 76, 210–236. [Google Scholar] [CrossRef]
  50. Son, S.; Chen, L.; Kang, Q.; Derome, D.; Carmeliet, J. Contact Angle Effects on Pore and Corner Arc Menisci in Polygonal Capillary Tubes Studied with the Pseudopotential Multiphase Lattice Boltzmann Model. Computation 2016, 4, 12. [Google Scholar] [CrossRef] [Green Version]
  51. Zheng, J.; Ju, Y.; Wang, M. Pore-Scale Modeling of Spontaneous Imbibition Behavior in a Complex Shale Porous Structure by Pseudopotential Lattice Boltzmann Method. J. Geophys. Res. Solid Earth 2018, 123, 9586–9600. [Google Scholar] [CrossRef]
  52. Liu, H.; Kang, Q.; Leonardi, C.R.; Schmieschek, S.; Narváez, A.; Jones, B.D.; Williams, J.R.; Valocchi, A.J.; Harting, J. Multiphase Lattice Boltzmann Simulations for Porous Media Applications: A Review. Comput. Geosci. 2016, 20, 777–805. [Google Scholar] [CrossRef] [Green Version]
  53. Akai, T.; Alhammadi, A.M.; Blunt, M.J.; Bijeljic, B. Modeling Oil Recovery in Mixed-Wet Rocks: Pore-Scale Comparison between Experiment and Simulation. Transp. Porous Media 2019, 127, 393–414. [Google Scholar] [CrossRef] [Green Version]
  54. Guo, R.; Dalton, L.E.; Fan, M.; McClure, J.; Zeng, L.; Crandall, D.; Chen, C. The Role of the Spatial Heterogeneity and Correlation Length of Surface Wettability on Two-Phase Flow in a CO2-Water-Rock System. Adv. Water Resour. 2020, 146, 103763. [Google Scholar] [CrossRef]
  55. Cooper, L.J.; Daly, K.R.; Hallett, P.D.; Naveed, M.; Koebernick, N.; Bengough, A.G.; George, T.S.; Roose, T. Fluid Flow in Porous Media Using Image-Based Modelling to Parametrize Richards’ Equation. Proc. R. Soc. Math. Phys. Eng. Sci. 2017, 473, 20170178. [Google Scholar] [CrossRef] [Green Version]
  56. Czachor, H. Modelling the Effect of Pore Structure and Wetting Angles on Capillary Rise in Soils Having Different Wettabilities. J. Hydrol. 2006, 328, 604–613. [Google Scholar] [CrossRef]
  57. Murison, J.; Semin, B.; Baret, J.-C.; Herminghaus, S.; Schröter, M.; Brinkmann, M. Wetting Heterogeneities in Porous Media Control Flow Dissipation. Phys. Rev. Appl. 2014, 2, 034002. [Google Scholar] [CrossRef] [Green Version]
  58. Benard, P.; Kroener, E.; Vontobel, P.; Kaestner, A.; Carminati, A. Water Percolation through the Root-Soil Interface. Adv. Water Resour. 2016, 95, 190–198. [Google Scholar] [CrossRef]
  59. Cassie, A.B.D. Contact Angles. Discuss. Faraday Soc. 1948, 3, 11–16. [Google Scholar] [CrossRef]
  60. Patel, R.A. Yantra. Available online: https://bitbucket.org/yantralbm/yantra (accessed on 1 January 2022).
  61. Patel, R.A. Lattice Boltzmann Method Based Framework for Simulating Physico-Chemical Processes in Heterogeneous Porous Media and Its Application to Cement Paste. Ph.D. Thesis, Ghent University, Ghent, Belgium, 2016. [Google Scholar]
  62. Sukop, M.C.; Or, D. Lattice Boltzmann Method for Modeling Liquid-Vapor Interface Configurations in Porous Media: Lattice Boltzmann for Interface Configurations. Water Resour. Res. 2004, 40, W01509. [Google Scholar] [CrossRef]
  63. Sukop, M.C.; Thorne, D.T. Lattice Boltzmann Modeling: An Introduction for Geoscientists and Engineers; Springer: Berlin/Heidelberg, Germany; New York, NY, USA, 2006; ISBN 978-3-540-27981-5. [Google Scholar]
  64. Wen, B.; Huang, B.; Qin, Z.; Wang, C.; Zhang, C. Contact Angle Measurement in Lattice Boltzmann Method. Comput. Math. Appl. 2018, 76, 1686–1698. [Google Scholar] [CrossRef] [Green Version]
  65. Liu, H.; Valocchi, A.J.; Kang, Q.; Werth, C. Pore-Scale Simulations of Gas Displacing Liquid in a Homogeneous Pore Network Using the Lattice Boltzmann Method. Transp. Porous Media 2013, 99, 555–580. [Google Scholar] [CrossRef]
  66. Diotallevi, F.; Biferale, L.; Chibbaro, S.; Lamura, A.; Pontrelli, G.; Sbragaglia, M.; Succi, S.; Toschi, F. Capillary Filling Using Lattice Boltzmann Equations: The Case of Multi-Phase Flows. Eur. Phys. J. Spec. Top. 2009, 166, 111–116. [Google Scholar] [CrossRef]
  67. Pooley, C.M.; Kusumaatmaja, H.; Yeomans, J.M. Modelling Capillary Filling Dynamics Using Lattice Boltzmann Simulations. Eur. Phys. J. Spec. Top. 2009, 171, 63–71. [Google Scholar] [CrossRef]
  68. Lu, G.; Wang, X.-D.; Duan, Y.-Y. Study on Initial Stage of Capillary Rise Dynamics. Colloids Surf. Physicochem. Eng. Asp. 2013, 433, 95–103. [Google Scholar] [CrossRef]
  69. Sun, C.; McClure, J.E.; Mostaghimi, P.; Herring, A.L.; Berg, S.; Armstrong, R.T. Probing Effective Wetting in Subsurface Systems. Geophys. Res. Lett. 2020, 47, e2019GL086151. [Google Scholar] [CrossRef]
  70. McHale, G.; Newton, M.I.; Shirtcliffe, N.J. Water-Repellent Soil and Its Relationship to Granularity, Surface Roughness and Hydrophobicity: A Materials Science View. Eur. J. Soil Sci. 2005, 56, 445–452. [Google Scholar] [CrossRef] [Green Version]
  71. Carmeliet, J.; Chen, L.; Kang, Q.; Derome, D. Beyond-Cassie Mode of Wetting and Local Contact Angles of Droplets on Checkboard-Patterned Surfaces. Langmuir 2017, 33, 6192–6200. [Google Scholar] [CrossRef]
  72. Zhang, B.; Wang, J.; Liu, Z.; Zhang, X. Beyond Cassie Equation: Local Structure of Heterogeneous Surfaces Determines the Contact Angles of Microdroplets. Sci. Rep. 2015, 4, 5822. [Google Scholar] [CrossRef] [Green Version]
  73. Liu, H.; Ju, Z.; Bachmann, J.; Horton, R.; Ren, T. Moisture-Dependent Wettability of Artificial Hydrophobic Soils and Its Relevance for Soil Water Desorption Curves. Soil Sci. Soc. Am. J. 2012, 76, 342–349. [Google Scholar] [CrossRef]
  74. Landl, M.; Phalempin, M.; Schlüter, S.; Vetterlein, D.; Vanderborght, J.; Kroener, E.; Schnepf, A. Modeling the Impact of Rhizosphere Bulk Density and Mucilage Gradients on Root Water Uptake. Front. Agron. 2021, 3, 622367. [Google Scholar] [CrossRef]
Figure 1. For the calibration of the contact angle the dependence on the pseudo density of the wall ρ w a l l was determined with the help of simulations. Illustrated are the exemplary results for CA 30°, 90° and 150° with the liquid (blue) and solid phases (grey).
Figure 1. For the calibration of the contact angle the dependence on the pseudo density of the wall ρ w a l l was determined with the help of simulations. Illustrated are the exemplary results for CA 30°, 90° and 150° with the liquid (blue) and solid phases (grey).
Water 14 01110 g001
Figure 2. Dependency of equilibrium CA on pseudo wall density ρ w a l l with fitting line and equation.
Figure 2. Dependency of equilibrium CA on pseudo wall density ρ w a l l with fitting line and equation.
Water 14 01110 g002
Figure 3. Simulation setup for capillary intrusion test of water (blue) into parallel plates.
Figure 3. Simulation setup for capillary intrusion test of water (blue) into parallel plates.
Water 14 01110 g003
Figure 4. Comparison of the time-dependent penetration length for LBM and analytical Washburn (WB) equation at CA 0° and 60°.
Figure 4. Comparison of the time-dependent penetration length for LBM and analytical Washburn (WB) equation at CA 0° and 60°.
Water 14 01110 g004
Figure 5. Comparison of the time-dependent penetration length for LBM and analytical Washburn equation with alternating CA coating for the CA pairs: 0° and 30° as well as 0° and 60° with an equidistant strip width of 0.2   mm .
Figure 5. Comparison of the time-dependent penetration length for LBM and analytical Washburn equation with alternating CA coating for the CA pairs: 0° and 30° as well as 0° and 60° with an equidistant strip width of 0.2   mm .
Water 14 01110 g005
Figure 6. Simplified porous medium consisting of circular particles (orange) and water source (blue) at the left side.
Figure 6. Simplified porous medium consisting of circular particles (orange) and water source (blue) at the left side.
Water 14 01110 g006
Figure 7. Coating patterns of the simplified porous medium: (a) homogeneous, (b) body-throat with a specific CA of each compartment type where the extension of the compartment with reduced wettability is defined by angle α and (c) stripes defined by angle α . Figure part (b) shows the throat-phobic coating with α = 30 ° and (c) the stripes coating with α = 15 ° .
Figure 7. Coating patterns of the simplified porous medium: (a) homogeneous, (b) body-throat with a specific CA of each compartment type where the extension of the compartment with reduced wettability is defined by angle α and (c) stripes defined by angle α . Figure part (b) shows the throat-phobic coating with α = 30 ° and (c) the stripes coating with α = 15 ° .
Water 14 01110 g007
Figure 8. Simplified three-dimensional porous medium, consisting of spherical particles (brown) with reduced wettable coating (red) and water source (blue) at the left side.
Figure 8. Simplified three-dimensional porous medium, consisting of spherical particles (brown) with reduced wettable coating (red) and water source (blue) at the left side.
Water 14 01110 g008
Figure 9. Sandy soil with mucilage coating (with mucilage: CA = 100°, without: CA = 20°).
Figure 9. Sandy soil with mucilage coating (with mucilage: CA = 100°, without: CA = 20°).
Water 14 01110 g009
Figure 10. Infiltration for homogeneous coating with CA = 50 ° at (a) simulation start, (b) 0.5   ms , (c) 3   ms as well as (d) the equilibrium state for homogeneous coating with CCA (red) 54° for > 0.5   ms .
Figure 10. Infiltration for homogeneous coating with CA = 50 ° at (a) simulation start, (b) 0.5   ms , (c) 3   ms as well as (d) the equilibrium state for homogeneous coating with CCA (red) 54° for > 0.5   ms .
Water 14 01110 g010
Figure 11. Penetration length as function of time for different CA at homogeneous coating.
Figure 11. Penetration length as function of time for different CA at homogeneous coating.
Water 14 01110 g011
Figure 12. v ¯ as function of CA for homogenous coating simulated by LBM.
Figure 12. v ¯ as function of CA for homogenous coating simulated by LBM.
Water 14 01110 g012
Figure 13. v ¯ as a function of CA for body-throat coating pattern: (a) body-phobic and (b) throat-phobic. The CA is related to the compartments of the particle surface whose wettability is reduced.
Figure 13. v ¯ as a function of CA for body-throat coating pattern: (a) body-phobic and (b) throat-phobic. The CA is related to the compartments of the particle surface whose wettability is reduced.
Water 14 01110 g013
Figure 14. Average infiltration velocity v ¯ as a function of CA for body-throat coating patterns to compare the effects of body-phobic and throat-phobic coatings. The CA refers to those compartments of the particle surface whose wettability is reduced.
Figure 14. Average infiltration velocity v ¯ as a function of CA for body-throat coating patterns to compare the effects of body-phobic and throat-phobic coatings. The CA refers to those compartments of the particle surface whose wettability is reduced.
Water 14 01110 g014
Figure 15. v ¯ as a function of CA for stripes coating. The CA is related to the compartments of the particle surface whose wettability is reduced.
Figure 15. v ¯ as a function of CA for stripes coating. The CA is related to the compartments of the particle surface whose wettability is reduced.
Water 14 01110 g015
Figure 16. v ¯ as a function of CA for homogeneous, body-phobic, throat-phobic and stripes coating pattern. The CA is related to the compartments of the particle surface whose wettability is reduced.
Figure 16. v ¯ as a function of CA for homogeneous, body-phobic, throat-phobic and stripes coating pattern. The CA is related to the compartments of the particle surface whose wettability is reduced.
Water 14 01110 g016
Figure 17. Infiltration for body-phobic coating with CA = 70 ° at (a) 0.5   ms , (b) 2.5   ms and for throat-phobic coating with CA = 70 ° at (c) 0.5   ms , (d) 2.5   ms as well as homogeneous coating with corresponding effective Cassie CA = 48 ° at (e) 0.5   ms and (f) 2.5   ms . The infiltrating water is shown in blue, while the particle surface with CA = 0 ° is orange, with CA = 70 ° is green and with effective Cassie CA = 48 ° is brown.
Figure 17. Infiltration for body-phobic coating with CA = 70 ° at (a) 0.5   ms , (b) 2.5   ms and for throat-phobic coating with CA = 70 ° at (c) 0.5   ms , (d) 2.5   ms as well as homogeneous coating with corresponding effective Cassie CA = 48 ° at (e) 0.5   ms and (f) 2.5   ms . The infiltrating water is shown in blue, while the particle surface with CA = 0 ° is orange, with CA = 70 ° is green and with effective Cassie CA = 48 ° is brown.
Water 14 01110 g017
Figure 18. Comparison of coatings with 1:1 ratio of surface areas with hydrophilic and reduced wettability as well as the corresponding homogeneous coating according to the Cassie equation.
Figure 18. Comparison of coatings with 1:1 ratio of surface areas with hydrophilic and reduced wettability as well as the corresponding homogeneous coating according to the Cassie equation.
Water 14 01110 g018
Figure 19. CCA for all coatings as a function of α .
Figure 19. CCA for all coatings as a function of α .
Water 14 01110 g019
Figure 20. Infiltration in three-dimensional water conductive system with surface stripe width of α = 10 ° at (a) 0.1   ms , (b) 0.2   ms and (c) 0.7   ms as well as (d) equilibrium state of the water repellent system with α = 15 ° reached at 0.1   ms .
Figure 20. Infiltration in three-dimensional water conductive system with surface stripe width of α = 10 ° at (a) 0.1   ms , (b) 0.2   ms and (c) 0.7   ms as well as (d) equilibrium state of the water repellent system with α = 15 ° reached at 0.1   ms .
Water 14 01110 g020
Figure 21. Soil system with (a) homogeneous CA = 20 is water conductive, while (b) with different CA (with mucilage: CA = 100 ° , without: CA = 20 ° ) is water repellent. Both images show infiltration after 11 ms, while in (b) the equilibrium state is present, which is already reached after 0.5 ms.
Figure 21. Soil system with (a) homogeneous CA = 20 is water conductive, while (b) with different CA (with mucilage: CA = 100 ° , without: CA = 20 ° ) is water repellent. Both images show infiltration after 11 ms, while in (b) the equilibrium state is present, which is already reached after 0.5 ms.
Water 14 01110 g021
Table 1. Listing and description of all basic particle coating patterns used for the simplified porous medium.
Table 1. Listing and description of all basic particle coating patterns used for the simplified porous medium.
Coating PatternDescription
homogeneous patternall particles homogeneously coated with the same CA (Figure 7a)
body-throat patternparticle surface of the pore throat and pore body having different CA (Figure 7b)
stripes patternparticle surface divided into stripes of equal size, which alternately have two different CA (Figure 7c)
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Bentz, J.; Patel, R.A.; Benard, P.; Lieu, A.; Haupenthal, A.; Kroener, E. How Heterogeneous Pore Scale Distributions of Wettability Affect Infiltration into Porous Media. Water 2022, 14, 1110. https://doi.org/10.3390/w14071110

AMA Style

Bentz J, Patel RA, Benard P, Lieu A, Haupenthal A, Kroener E. How Heterogeneous Pore Scale Distributions of Wettability Affect Infiltration into Porous Media. Water. 2022; 14(7):1110. https://doi.org/10.3390/w14071110

Chicago/Turabian Style

Bentz, Jonas, Ravi A. Patel, Pascal Benard, Alice Lieu, Adrian Haupenthal, and Eva Kroener. 2022. "How Heterogeneous Pore Scale Distributions of Wettability Affect Infiltration into Porous Media" Water 14, no. 7: 1110. https://doi.org/10.3390/w14071110

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop