How Heterogeneous Pore Scale Distributions of Wettability Affect Inﬁltration into Porous Media

: Wettability is an important parameter that signiﬁcantly determines hydrology in porous media, and it especially controls the ﬂow of water across the rhizosphere—the soil-plant interface. However, the inﬂuence of spatially heterogeneous distributions on the soil particles surfaces is scarcely known. Therefore, this study investigates the inﬂuence of spatially heterogeneous wettability distributions on inﬁltration into porous media. For this purpose, we utilize a two-phase ﬂow model based on Lattice-Boltzmann to numerically simulate the inﬁltration in porous media with a simpliﬁed 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 speciﬁc 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 speciﬁc microscale locations of exudate depositions may improve models of root water uptake.


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 nonwettable (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.

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.

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 ), (1) 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.

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 where f i represents the particle distribution function along direction i in the velocity space, ∆t is the time step, τ NS 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 eq i is the equilibrium distribution function given by 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 and the macroscopic kinematic viscosity (ν) is related to the relaxation time as 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 int 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 int is given as where G = −120 is the interaction strength. ψ is the interaction potential given as per the Shan-Chen model [47]: 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 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 ρ wall ( 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.

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]: 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: 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.   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.

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): Table 1. Listing and description of all basic particle coating patterns used for the simplified porous medium.

Coating Pattern Description
homogeneous pattern all particles homogeneously coated with the same CA (Figure 7a • 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 yand 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 • . 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.

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 synchrotronbased X-ray tomographic microscopy (SRXTM). Mucilage structures were visible as twodimensional 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.

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 (Figures 10 and 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 • (Figures 10 and 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].

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). Figure 13. v as a function of CA for body-throat coating pattern: (a) body-phobic and (b) throatphobic. The CA is related to the compartments of the particle surface whose wettability is reduced.
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 α body-phob ≥ α throat-phob the effective total infiltration rate σ e f f corresponds to the infiltration rate of a system with only body-phobic coating. For α body-phob ≥ 20 • , σ e f f is even determined independently of α body-throat 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. 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.
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 threedimensional 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.

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). 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.
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. 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.
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.

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 (Figures 17 and 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 CAs. 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.

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.

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.

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.

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: With the help of the conversion factors of the dynamic viscosity C µ of fluid and the surface tension C γ the conversion factors for time and mass can be calculated by reshaping and inserting Equations (A2) and (A3): The kinetic viscosity in lattice units ν lb is determined by Equation (5) with τ NS = 1, which leads to ν lb = 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 µ lb,v = 14.31 mu lu ts and that of the fluid µ lb,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 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 γ.

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 with x and x f inest 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.