Evaporation Study Based on Micromodel Experiments: Comparison of Theory and Experiment

Evaporation—a key process for water exchange between soil and atmosphere—is controlled by internal water fluxes and surface vapor fluxes. Recent studies demonstrated that the dynamics of the water flow in corners determine the time behavior of the evaporation rate. The internal water flux of the porous media is often described by capillary flow assuming complete wetting. Particularly, the crucial influence of partial wetting, that is, the nonlinear contact angle dependency of the capillary flow has been neglected so far. The focus of the paper is to demonstrate that SiO2‐surfaces can exhibit contact angles of about 40°. This reduces the internal capillary flow by 1 order of magnitude compared to complete wetting. First, we derived the contact angle by inverse modeling. We conducted a series of evaporation experiments in a 2‐D square lattice microstructure connected by lognormal distributed throats. We used an explicit analytical power series solution of the single square capillary model. A contact angle of 38° ± 1° was derived. Second, we directly measured the contact angle of the Si‐SiO2 wafer using the Drop Shape Analyzer Krüss 100 and obtained an averaged contact angle of 42° ± 2°. The results support the single square capillary model as an appropriate model for the description of the evaporation process in an ideal square capillary.


Introduction
Evaporation is a key process for water exchange between soil and atmosphere. This complex multiphase process is determined by the diffusive transport through the soil gas phase and the mass transfer boundary layer as well as the liquid water transport through the porous media. The evaporation rate exhibits a two-stage behavior that is often experimentally observed under low atmospheric demand (<5 mm/day; Shahraeeni et al., 2012). Stage 1 is characterized by a constant evaporation rate and stage 2 by a falling evaporation rate (Brutsaert & Chen, 1995). During stage 1, enough water (driven by the atmospheric demand/atmospheric conditions) is transported from the drying front (interface between saturated and unsaturated/drying soil) to the soil surface via a continuous water phase cluster, that is, through pores which are hydraulically connected.
The main water transport mechanism is capillary flow (capillary pumping), where the pore structure and microstructure of the solid surface (inner surface roughness) determine the water flow from the pore corners (corner flow; Figure 1) or from the surface cavities (dimension of about 10 μm; thick-film flow) to the drying front (Blunt & Scher, 1995;Hashemi et al., 1999;Lehmann et al., 2008;Sahimi, 2011). The complex interplay of both flow types controls the kinetic interphase mass transfer near the soil surface, where the driving vapor concentration gradient is determined by atmospheric conditions (wind velocity v ∞ , temperature T, relative humidity RH = f (p H2O ), p H2O -water vapor pressure), which are usually expressed by one lumped parameter-the viscous boundary layer thickness δ = f(v ∞ , T, RH) (Figure 1; Schlichting & Gersten, 2017;Geistlinger & Leuther, 2018;Peters et al., 2015;Lehmann et al., 2008;Yiotis et al., 2007;Prat, 2002).
Water flows in (i) pore channels (denoted as duct or bulk flow), (ii) corners (corner flow), (iii) along surface cavities/capillaries (thick-film flow), and (iv) within thin water films at smooth pore walls (thin-film flow). The hierarchy of the flow sequence is determined by the water pressure gradients, given by the difference of water pressure of the bulk phase and the negative capillary pressure at the flow front (Lenormand & Zarcone, 1984). For a water-wet porous medium thin-film flow starts first (i.e., a SiO 2 surface is always covered by a thin water film; Kibbey, 2013) followed by thick-film flow, corner flow, and bulk flow (Geistlinger & Mohamadian, 2015). It has been shown that the contribution of thin-film flow to evaporation from rectangular capillaries is negligible compared to diffusion and corner flow (Eijkel et al., 2005). Note that the term thick-film flow is often used to describe corner flow (Lehmann et al., 2008;Yiotis et al., 2012a).

Modeling
There are two groups of evaporation models (i) pore network models (Prat, 2002(Prat, , 2011 and (ii) semiphenomenological models which couples the pore scale to the representative elementary volume scale under certain assumptions. The focus of this paper is on semiphenomenological models. The advantage of these models is that macroscopic, effective properties of the porous media can be measured (retention curve, unsaturated hydraulic conductivity, and evaporation rate) using the multistep outflow method (Durner & Iden, 2011;Vogel et al., 2008) and the HYPROP method (Schindler et al., 2010). This allows to study the complex interplay and the impact of certain pore scale properties on the interphase mass transfer near the surface, that is, on the evaporation rate.
Two important semiphenomenological models are shown in Figure 1: (a) the capillary bundle model of evaporation (isolated pore evaporation [IPE] model) that is based on the pore size distribution (PSD) of the porous medium and (b) the single square capillary model (SSC model) of evaporation that is based on corner flow driven by the high capillary pressure of water-filled corners inside a rectangular geometry of a waterwet pore throat or pore body. (Figure 1a). The IPE model couples the 2-D surface model of evaporation (Suzuki & Maeda, 1968) to the internal water flux neglecting viscous forces. The 2-D surface model of evaporation captures two key features of an evaporating water-saturated porous surface: (1) the discrete nature of the vaporization surface consisting of numerous isolated pores linked by capillarity with internal liquid reservoir (therefore IPE model) and (2) the nonlinear vapor diffusive fluxes from isolated pores across air viscous sublayer to the atmosphere. The mass transfer through the diffusive boundary layer, composed the individual mass fluxes of the isolated pores, is highly nonlinear. Interference between neighboring active pores results in a reduction of the diffusive flux from each pore compared to isolated pores, whereas the total flux from the entire Figure 1. Semiphenomenological models of evaporation: (a) the capillary bundle model is based on capillary pumping (q w ) between active pores represented by the small r 1 capillary and inactive pores represented by the large r 2 capillary. Evaporation lengths (drying front depths) L and L c are shown at two different times t and t c , where L c is the characteristic length just before the meniscus recedes into the small capillary (= detachment time; K(θ)-hydraulic conductivity of the porous medium represented by the orange block, e(t)-evaporation rate, δ-thickness of the viscous boundary layer, capillary heights h c1 and h c2 , gravitational length L G = h c1h c2 , and PSD-pore size distribution). (b) The single square capillary model is based on the steady state between corner flow q c within a square capillary (width d 0 ) and diffusion flow q diff (= evaporation rate density). The capillary is initially saturated and after a critical time t c (length of stage 1) the tips of the corner flow detach from the surface and recede into the capillary. At this critical time the "spanning" water cluster which connects the surface with duct flow region detaches from the surface. The x-extension of the spanning water cluster is given by the averaged depth of the percolation front 〈x p 〉(t c ) shown in Figure 1b. The corner flow region (sometime denoted as film region) extends from the percolation front x p (y) (duct flow front) and the evaporation front x i (y) (corner flow front). surface increases, since diffusion in lateral direction is suppressed by the "communicating" boundary condition (Lord Rayleigh, 1945). During evaporation the density of active pores decreases, which weakens lateral interference between pores. The diffusive flux of an individual active pore increases, whereas the total flux from the entire surface decreases. This nonlinearity of the evaporation rate is due to the dependence of the mass transfer rate constant on the ratio of boundary layer thickness and active pore size. The internal water flux of the porous media is realized by hydraulic connected region of the PSD. The PSD can be derived from the retention curve (capillary pressure-saturation relationship; Lehmann et al., 2008) or measured directly by X-ray micro-tomography and image analyses (Geistlinger & Leuther, 2018). The IPE model is based on the hypothesis that the hydraulically active 3-D pore network can be approximated by a bundle of active and inactive capillaries and that its size distribution can be described by the PSD of the porous medium. Based on the coupling of the 2-D surface model of evaporation (Suzuki & Maeda, 1968) and the internal water flux, the IPE model describes the whole stage 1 period of evaporation. We previously tested the IPE model for real soils using a complete set of hydraulic functions measured by the HYPROP method and found reasonable agreement (Geistlinger & Leuther, 2018). The critical point was to estimate the hydraulic connected region of the PSD.

IPE model
SSC model (Figure 1b). A limitation of the IPE model is that it is based on the curvature of the capillary menisci, that is, only duct flow is considered and corner flow is neglected. Previous work by Lenormand and Zarcone (1984) solved the corner flow problem using the hydraulic diameter approximation. Dong and Chatzis (1995) solved the corner flow problem for a general rectangular geometry (corner angle: 0 < α < 90°) and for partial wetting behavior (contact angle: 0 < θ < 90°). Based on this work, Yiotis et al. (2012a) coupled corner flow to gaseous diffusion flow and solved the quasi steady state problem (SSC model), that is, corner flow q c equalize diffusion flow q diff at the mean evaporation front 〈x i 〉 (Figure 1b). Note that the SSC model is based on the mean pore size defined as mean capillary width d 0 (Figure 1b), whereas the IPE model is based on the whole PSD (insert in Figure 1a). Starting with a complete water-filled capillary, evaporation needs some time to empty the capillary. The time needed until the water front (hydraulically connected cluster) detaches from the surface determines the detachment time t c and the length of the stage 1 period. The SSC model describes the evaporation rate also during the stage 2 period, that is, gaseous diffusion flow from the receding mean evaporation front.
Hypotheses of the SSC model. The quasi steady state between corner flow and gaseous diffusion flow determines evaporation during stage 1 and stage 2. The mean pore size defines the effective capillary of the porous media, and the detachment time determines the length of stage 1 evaporation.
An experimental test of the SSC model was conducted by Yiotis et al. (2012b) using glass beads packs. However, the test remained inconclusive in several points: 1. Rectangular geometry. That is, how to map the highly irregular shape of pore channels, which consist of convex-and concave-curved surface areas to a rectangular geometry of a square capillary. 2. Complete wetting. The theoretical treatment assumes unrealistic complete wetting (zero contact angle) for the glass beads (the fluid was pentane and not water; the contact angle of pentane on glass is about 30°; Li et al., 2015). It has been shown that contact angles on glass surfaces using standard cleaning methods lie between 15°and 30° (Iglauer et al., 2014) and that heterogeneous wetting properties can have a significant influence on capillary trapping (Geistlinger & Ataei-Dadavi, 2015). 3. Roundness factor. It is not clear why the irregularly shaped pore channels, which are often approximated by triangular cross sections (Blunt, 2001), should exhibit round corners. To obtain reasonable agreement between experiment and theory, Yiotis et al. (2012b) used a "roundness factor" as a fitting factor. However, it should be noted that a finite roundness factor will weaken the capillary pressure of the corners and has therefore the same influence as a finite contact angle. For instance, roundness factors of 0.4 and 0.7 (definition after Dong & Chatzis, 1995) correspond to contact angles of sharp corners of 12°and 63°, where the cross section of the wetting phase inside the round corner and sharp corner is equal.
The main objective of this paper is to test the fundamental hypothesis of the SSC model by micromodel experiments using Si-microfabricated micromodels with lognormal-distributed "ideal" square channels. An interval-based inductively coupled plasma-deep reactive ion etching (ICP-DRIE) technology (Küchler et al., 2003) for microfabrication was applied that results in high edge steepness and a true mapping of the lattice structure with depth (2°deviation from the vertical line).
A second objective is to derive the SSC model consequently as a function of the contact angle (θ) and temperature (T) to study the contact angle and temperature dependence on the evaporation process. Particularly, the study focuses on the impact of partial wetting on the evaporation flux at the surface of a porous structure, for example, the pore network of a soil column.
A section of the square lattice is shown in Figure 2b. Lattice sites are also called pores, and the half throat width is denoted as throat radius according to Yiotis et al. (2012a). The open side of the microstructure (red line in Figure 2a) corresponds to the top side in the figures of the paper.
A microfabrication method involving photolithography, ICP-DRIE, and anodic bonding was used to fabricate the micromodel in a silicon wafer. An interval-based ICP-DRIE technology (Küchler et al., 2003;Zue et al., 2013) was applied for anisotropic etching. This resulted in high edge steepness and a true mapping of the lattice structure with depth ( Figure 2b), that is, only minimal underetching is produced (1.2°deviation from the vertical line and maximal deviation from the horizontal bottom line of 10 μm). This 1.2°deviation was achieved for a ninefold deeper microstructure compared to the microstructure produced by Zue et al. (2013). Standard methods, like isotropic plasma etching and chemical etching (Vorhauer et al., 2015), lead to underetching and corner rounding effects, which have to be taken into account for the calculation of the capillary corner flow (roundness factor; Dong & Chatzis, 1995;Yiotis et al., 2012a).
Briefly, the microfabrication sequence is described as follows: a chromium-coated soda lime photomask of 7 in. size was patterned based on the Matlab-generated data. Next, a cleaned p-type silicon wafer with (100) orientation, 150-mm diameter, 675 μm thick, was covered with a 2-μm-thick oxide film by wet thermal oxidation and subsequently coated with approximately 8-μm-thick resist (AZ9260). The mask pattern was transferred from the photomask onto the resist using contact lithography and then transferred into the oxide layer by plasma etching on an SPTS DSi Rapier ICP etcher.
The patterned film stack served as a hard mask for the subsequent deep anisotropic silicon etch. The silicon trenches were etched to a nominal depth of 300 μm by DRIE carried out on the same DSi Rapier ICP etch tool. After DRIE, remaining resist was removed by oxygen plasma ashing and the oxide layer by buffered oxide etching. The structured silicon wafer was cleaned and thermally oxidized by dry oxidation to form a 50-nm-thin silicon oxide layer onto the silicon surface to ensure a similar wetting behavior like siliceous surfaces. The structured wafer was then covered by a borosilicate glass wafer (Schott Borofloat 33, 150-mm diameter, 1 mm thick) by anodic bonding (400°C, 600-800 V). One side of the glass wafer was cut such that it fitted the open side of the microstructure. This wafer stack was sawn to fit into a 144 × 100-mm 2 rectangle.

Water Resources Research
The schematic of the micromodel is shown in Figure 2a and the rectangular corner shape of the lattice structure in Figure 2b.

Contact Angle Measurement
The contact angles were measured at the left and right sides of a droplet of deionized water both for the Si-SiO 2 -wafer and for the Borofloat cover glass using the Drop Shape Analyzer Krüss DSA 100. The Si-SiO 2 -wafer was treated as the micromodel, that is, (i) ICP-DRIE-etching of the Siwafer and (ii) deposition of a 50-nm thermal SiO 2 layer. The experimental values are listed in Table 1. The mean contact angle was 48 ± 2°(Si-SiO 2wafer) and is 23 ± 2°(Borofloat-cover glass).

Methodology and Experimental Setup
We conducted three vertical evaporation experiments at three different temperatures: 28.3°C (Experiment 1), 41.5°C (Experiment 2), and 60.8°C (Experiment 3) and one horizontal experiment (Experiment 4) at 20.3°C. The ambient boundary conditions and the experimental results are listed in Table 2. Prior to each experiment, the micromodel was saturated with deionized water under vacuum for 2 hr using an evacuator. During the evaporation/drying process, the micromodel was positioned vertically (open boundary on top) or horizontally, and the mass of the micromodel and heating plate was measured by a high-precision digital balance (Sartorius Secura 1103-1S, x ± 0.001g) and recorded every 10 s by DasyLab software (version DasyLab 2016).
The micromodel was mounted on a heating metal plate (meandering heaters at the backside), which ensures a homogeneous temperature distribution over the micromodel with a spatial deviation (midpoint to boundary point) of 0.5°C for Exp. 1 and Exp. 4, 1.6°C for Exp. 2, and 4.2°C for Exp. 3. The temperature control works with a 1°tolerance. The surface temperature was measured with an infrared beam device (Trotec BP25). The relative humidity in the lab during experiments was recorded every 5-10 min by a data logger (dataTaker DT80 series2). The temperature difference between the micromodel and the heating plate was taken into account. Images of the fluid distribution were monitored by a high-resolution digital single-lens reflex camera (Canon 5-D Mark IV, lens: Canon EF 100mm F2.8L Marco IS USM) with a resolution of 6,720 × 4,480 pixels at process-dependent time intervals (early stage 1: 1 min; late stage 2: 30 min). To enhance the contrast of different fluidic phases, a circular cold lighting source was applied in front of the micromodel. A miniphotographic studio (white walls around the micromodel; see Figure 3) was used to produce homogeneous light conditions and to minimize air convection in the lab.

Image Processing
The PSD and the fluid patterns were analyzed using Fiji/ImageJ (Schindelin et al., 2012).

PSD
The PSD of the micromodel was analyzed using the maximum-inscribed sphere method implemented in the ImageJ plugin "Local thickness." The PSD is shown in Figure 4 (blue dots). The best fit curve (lognormal distribution, red curve) yields a mean pore radius of 197 μm with a standard deviation of 55 μm. Note  Note. Notation: Bo-bond number, T-temperature in degrees Celsius, (RH)-relative humidity, θ (t c )-contact angle derived from t c , θ (x p )-contact angle derived from x p , j diff , IP-evaporation flux density derived from image processing (IP), δ 0 -boundary layer thickness defined by the molecular vapor diffusion coefficient.
that the throat distribution is characterized by a mean throat radius of 200 μm and a standard deviation of 31 μm. Since the pores located at the nodes (lattice sites) are always larger than the four neighboring throats, the PSD exhibits a broader maximum and is shifted to larger pore sizes causing a higher standard deviation.

Fluid Pattern Analysis
1. Image registration. Prior to image acquisition during the evaporation process, a grains mask image was taken in dry state to conduct some image calculation. The registration plugin "TurboReg" was applied to align the grain image with the fluid pattern image. 2. Brightness normalization. The fluid pattern image exhibited a radial illumination drift because of the circular lighting source. To normalize the brightness, the mean gray value was subtracted as a function of radius from the image using the "Math" function yielding a corrected image with homogenous illumination. 3. Phase separation. Lights irradiating at grains edges reflected back along its original way in dry pores and water saturated pores. Thus, bright lines along grain walls in dry pores were observed in the raw RGB image. In film pores the slope of corner flow eliminates this light artifact. After the brightness was normalized, a binary image was computed by the autothreshold method "Otsu." Subsequently, wall artifacts were removed by applying the "Dilate"-"Fill Holes" operations twice followed by two times the "Erode" operation. The "Median" filter was applied followed by an "Area opening" with the MorphoLibJ plugin (Legland et al., 2016) to eliminate segmentation noise. 4. Fluid pattern analysis. Water saturation was obtained by the "Analyze Particles" tool summarizing the total pixels of water clusters. The functions "getDimensions" and "getStatistics" yielded the coordinates of each pixel on the percolation front/evaporation front.  . Pore size distribution (blue dots) of the micromodel and best fit (red curve, lognormal distribution). The black curve shows the throat size distribution. The inset shows the spatial distribution of pore sizes for a 4 × 4 lattice section (orange: large radius, purple: small). PSD = pore size distribution.

Theoretical Background: 1-D Analytical Theory of the SSC Model
In this section we derive analytical expressions for the time-dependent mass loss, the evaporation front, and the percolation front. This will be used for the comparison to the experimental results in the next section.
Evaporation from a saturated porous media is controlled by capillary flows from large to small capillaries. Hence, the corner flow and water menisci Σ across the pore channel occur first in large capillaries. The basic assumptions of the SSC model are that (i) the gas-water interface near the meniscus Σ is stationary and (ii) that the curvature along the x direction can be neglected when compared to the curvature in y direction (rectangular to the x direction). Hence, the capillary pressure p c of the water meniscus at x = −x Σ must be equal to the capillary pressure of the gas-water interface at x = −x p (Figure 5a) where H denotes the curvature and r x and r y the radii in principal directions x and y, respectively, θ the contact angle, and σ gw (N/m) the interfacial tension of the gas-water interface. The radius r y is (Dong & Chatzis, 1995;Legait, 1983) where r 0 is half of the mean capillary width (in the following denoted as radius). Equation (2a) determines the validity range of the SSC model, that is, the following discussion considers only the range x ≥ −x p . The θ dependence of the corner flow Q wc (m 3 /s; shown by the blue arrows in Figure 5a) is determined by two factors: (i) a reduced cross section A θ ð Þ ¼ C * θ ð Þr 2 y , (m 2 ) shown in Figure 5b and (ii) a higher flow resistance β(θ) (-) with the integral

10.1029/2018WR024647
Water Resources Research μ w (Pa·s) denotes the dynamic viscosity of the water phase, α θ ð Þ ¼ β θ ð Þ C * θ ð Þ the rescaled flow resistance (see Figure 6), ρ = r y /r p the normalized radius of curvature, and Bo ¼ ρ w gr 2 p =σ gw the bond number. Typical bond numbers for the evaporation experiments are in the order of 10 −3 . The θ-dependent constant C* is given by (Dong & Chatzis, 1995) As the evaporation process is a slow drainage process (capillary number Ca = μ w u w /σ gw ≅ 4 × 10 −8 for Exp. 3, u w -velocity of the receding percolation front), the steady state between water flow and averaged evaporation flux 〈Q ev 〉 is always established. The mass balance over a small volume dV = A · dx (Figure 5a) at the evaporation front where m w denotes the water mass of the liquid phase in the volume element dV and ρ w the water density. The averaged evaporation flux is given by first Fick's law where A g denotes the cross section of the gas phase at the evaporation front. We approximate A g by A 0 g ¼ d 2 0 , because the water saturation near the surface (only four small tips of the corner flow) is approximately 0.
j diff x (kg/m 2 /s) denotes the mass evaporation rate density and D eff g (m 2 /s) the effective gas phase diffusion coefficient of the water molecules inside the porous media. For D eff g the Millington-Quirk tortuosity model (Millington & Quirk, 1961) is used, that is, D eff g ¼ D 0 g τ S g À Á (D 0 g -gas phase diffusion coefficient of the water molecules). The tortuosity τ = (ϕ) 1/3 (1 − S g ) 7/3 (ϕ-porosity, S g -gas saturation). The temperature depen- 1:81 m 2 =s (T in K; Shahraeeni et al., 2012). ΔC(x) = C 3 -C 2 is the difference of the vapor concentrations (kg/m 3 ) x 3 = δ (boundary layer thickness) and x 2 (mean evaporation front). The positive x direction is parallel to the normal surface vector ( Figure 5). The steady-state evaporation rate for stage 1 and stage 2 is determined by the following concentration gradients: where C a denotes the ambient vapor concentration at viscous boundary layer thickness δ and C s the saturated vapor concentration (kg/m 3 ). Note that the mass evaporation rate density j diff x (kg/m 2 /s) is related to volumetric evaporation rate density _ e by the factor ρ w (often denoted as evaporation rate, see Lehmann et al., 2008).
Inserting equations (3a) and (6a) into equation (5) yields where the new variable ϕ(x) is defined as Figure 6. Contact angle dependence of the corner flow caused by a reduced cross section and an increased flow resistance (inset). The red curve presents the best fit to the data points (blue dots). Those points are taken from Ransohoff and Radke (1988, Table 4) for the case of a square capillary (corner angle = 90°) and a free water surface (reduced viscosity = 0). The combined effect is expressed by rescaled flow resistance) α(θ) = β(θ)/C * (θ).

10.1029/2018WR024647
Water Resources Research The θand T-dependent constant b (right side equation (7a)) is the key variable of the SSC model where ν w denotes kinematic viscosity (m 2 /s) and T the temperature. Note that the kinematic viscosity and the surface tension exhibit a strong T dependence in the temperature range of interest, that is, between 20 and 60°C. Inside the capillary (inside the porous medium), b(θ) is zero (ΔC = 0) and ϕ(x) is determined by the Laplace equation The constant b is determined by the flux continuity boundary condition at the evaporation front x 2 which yields b = b(θ,T).
Before we derive an explicit solution of the vertical case (Bo ≠ 0), it is instructive to discuss the simpler horizontal case (Bo = 0). The steady-state solution ϕ(x) (equation (8b)) and its relation to the evaporation rate (equation (9)) is schematically shown in Figure 7. The concentration gradient shown on the right determines the slope of ϕ(x), that is, the θand T-dependent constant b(θ,T) (equation (7c)). With increasing time the evaporation front x 2 (t c + Δt) recedes into the porous medium and causes a weaker concentration gradient. The linear solution is "shifted" by Δx 1 and causes a change of the percolation front Δx 1 . It is obvious that the slope of the linear solution b(θ,T) gives the inverse of x 2 (t)-x 1 (t) (equation (12a)). Note that a small change of the evaporation front Δx 2 leads to a large change of the percolation front Δx 1, that is, the ratio Δx 1 /Δx 2 > 1 and increases with time. Physically, the decreasing evaporation flux leads to a decreasing corner flow (flux condition; decreasing viscous forces) which is proportional to dϕ dx ; hence, the slope of ϕ decreases as shown by the blue line in Figure 7 (left diagram). The white area above the red line (ϕ(x) diagram) determines the evaporated mass Δm w at t c and the blue area the increase of Δm w with time step Δt.

Water Resources Research
Now we consider the vertical case (Bo ≠ 0). Introducing a new variable ξ(x) = ρ 3 (x), we reformulate solution ((8b)) into a nonlinear integral equation for ξ where ε ¼ 3Bo rp . The explicit solution of this integral equation is obtained by a power series expansion in εΔx (Δx = x 2 -x), and the nth approximation of ξ(x) is given by with a k = 4/3 + k. Note that the solution for the horizontal case is given by ξ (0) (x), that is, by setting ε = 0. The convergence of the power series solution (11) is proven for n = 3 (Bo ≈ 10 −3 ) in Appendix B.
We now calculate measurable quantities of the evaporation process.

Distance between the evaporation front and percolation front
Due to the boundary condition ϕ(x 1 ) = 1 (equation (7b)), the distance between the evaporation front and percolation front, x 2 -x 1 , for horizontal flow (Bo = 0) is given by the inverse of the constant b(θ,T) For the vertical case (Bo ≠ 0) the implicit equation has to be solved.
At t = t c (x 2 = 0), equation (12a) yields the x extension (width) of the spanning water cluster x p (t c ).

Detachment time t c
The detachment time is given by the ratio of the evaporated mass (= mass loss) Δm w (white area above the red line in ϕ(x) diagram; Figure 7) and the stage 1 evaporation rate (equations (6a) and (6b)) 3. Time-dependence of the mass loss curve Δm w (t) during stage 2 We next derive the time dependence of the mass loss Δm w (t) and the evaporation front x 2 (t). We first discuss the analytical solution of the horizontal case and then describe the solution method for the vertical case. The calculation procedure is illustrated in Figure 7. To obtain the mass loss, we calculate the residual water mass inside the four corners m w (area below the blue curve in the ϕ(x)-diagram) and subtract this value from the water mass of the water-filled capillary m 0 (area of the dashed rectangle in the ϕ(x) diagram). The general expression for the mass loss is For the horizontal case we obtain a linear Δm w (x 2 ) relationship using solution (8b) and equation (6c) Δm w x 2 ð Þ ¼ a 3 −a 4 x 2 (14b) where the θand T-dependent constants a i are given in Appendix A. Given the time dependence x 2 (t), equation (14b) yields the time dependence of the mass loss. To derive the differential equation for x 2 (t), we consider the x 2 derivative of Δm w Inserting the x 2 dependence of the evaporation rate (equations (6a) and (6c)) yields the ordinary differential equation for x 2 (t) and its solution is The time dependence of the percolation front is obtained by inserting x 2 (t) into equation (12a) We discuss different limit cases of the steady-state solution. For 2 a6δ 2 ·Δt≪1 Δm w (t) shows a linear increase that extrapolates the linear time dependence of stage 1 for a small transition region between stage 1 and stage 2. For t >> 2 a6δ 2 ·Δt, Δm w (t) shows the expected ffiffi t p dependence.
To the best of our knowledge, we derived for the first time these explicit analytical expressions for the time dependence of the mass loss (equation (14a)), the evaporation front (equation (17a)), and the percolation front (equation (18)).
For the vertical case we follow the same solution procedure using the capability of Mathematica to define analytical functions f(x), where x is the solution of an implicit equation (e.g., equation (12b)) or the upper integration limit of a definite integral (e.g., equation (3b)).
It is instructive to compare the solutions of the horizontal and vertical evaporation processes (black and red curves in Figure 8). Given that the vertical flow length is shorter than the horizontal flow length, we expect that the evaporation front recedes faster into the porous media. We also expect that the extension (width) of the corner flow region is smaller compared to the horizontal case. This characteristic behavior is shown in Figures 8a and 8b. For the horizontal case more water can be transported to the evaporation front. Hence, the mass loss is higher compared to the vertical case as shown in Figure 8c.  Table 2). t c denotes the detachment time.

Water Resources Research
Note that the solutions satisfy the same boundary conditions, that is, yield the same width of the corner flow region at the critical time t c . (x p (t c ) = 4.2 mm; see  Figure 9a shows a typical spatial pattern of the corner flow region bounded by the evaporation (red) and percolation (blue) front during stage 2 at t = 244 min (Exp. 3; see Table 2). In Figure 9b (enlarged section of Figure 9a, yellow rectangle), the gas-water interface with the water menisci at the lattice nodes can be recognized (curved white lines). The yellow arrow shows a corner flow path from the percolation to the evaporation front, that is, through the unsaturated corner flow region. The shortest path cannot be taken, as many channels are blocked, that is, completely filled with water. The location of such a single-channel water cluster is indicated by the red arrow. These isolated single-channel or multichannel water cluster force the corner flow along a highly tortuous path. It is instructive to compare Figure 9b with Figure 10, which is an image-processed section of Figure 9b (yellow rectangle). The different phases are indicated by different colors: (i) the water phase of the duct flow area and the isolated water cluster are marked in blue, and (ii) the yellow areas indicate the vapor-saturated air demonstrating a continuous tortuous pathway of the vapor molecule through the corner flow region. The corner flow pathways are shown by the red areas. Most of the pathways are continuous pathways connecting the water-saturated region (blue region) with the completely dry region (yellow region), that is, the pathway i starts at a certain point (x i 1 , y i 1 ) at the heterogeneous percolation front (index 1) and ends at a certain point (x i 2 , y i 2 ) at the evaporation front (index 2).

Experimental Results and Discussion
Since the evaporation process is a drainage process, the pathway should always follow networks of connected channels with large width, that is, low capillary pressure. The water pressure gradient dissipates along such tortuous corner flow pathways. Besides continuous corner flow pathways, there are also isolated corner flow pathways shown in Figure 10 inside the dashed line square. Note that there is still a corner flow from the waterfilled channels (blue sections) to the tips of the corner flow pathway. Vorhauer et al. (2015) observed in regular stochastic Si-lattices that stable water rings around single grains during the evaporation process. Why are Figure 9. (a) Spatial stochastic pattern of the corner flow region bounded by the evaporation (red) and percolation (blue) front during stage2 at t = 244 min for evaporation Exp. 3 (see Table 2). (b) Magnified section of Figure 9a (yellow window). Water menisci, single isolated water-filled channels (red arrow), and two corner flow paths (yellow arrows) are shown. they absent in our experiments? For the following discussion we have to distinguish two types of corners: (i) corners along the flow direction and (ii) corner perpendicular to the flow denoted as grain corners. It is obvious that the water ring surface around the grain corner exhibits a saddle point geometry, that is, convex (positive) curvature within the corner and concave (negative) curvature around the grain. To keep the water pressure constant (p g -p w = constant = p c ), the negative curvature has to be compensated by a higher positive curvature. For sharp corners it is rather likely that these "tied-up" water rings become unstable and snap-off will occur. Therefore, isolated water-filled channels are more likely than water rings. The absence of water rings is an indirect indicator of the corner stiffness (etching precision) of the square lattice shown in Figure 2b. For round corners, the probability that water rings are thermodynamically stable is higher as discussed by Vorhauer et al. (2015).
Four different evaporation experiments with each about 20 snapshots of the fluid distribution are the basis for the quantitative analysis of the geometric characteristics of the evaporation process. We test the accuracy of the IP images by calculating the water saturation (summation over all water pixels) and by deriving the mass loss curve at about 20 time points. The relative error between the gravimetrical mass loss curve (thin black line in Figure 11a represents the best fit to the data) and the IP-derived mass loss (red dots in Figure 11) is 5% for Exp. 3 (11% for Exp. 1 and 9% for Exp. 2), that is, both data sets show good agreement.

Extension of the spanning water cluster x p (t c )
We start our discussion with stage 1 of the evaporation process. Figure 11a shows the mass loss curve for Exp. 3. The straight blue line indicates the stage 1 slope of the mass loss curve, which determines the constant evaporation rate 〈Q ev 〉 . This also determines the boundary layer thickness of 0.30 mm via equations (6a) and (6b) listed in Table 2. The experimental value of the extension (= width) of the spanning water cluster is 4.2 mm (Table 2). Solving the nonlinear equation (12a) for the contact angle (4.2 mm = 1/b(θ,T)) yields a contact angle of 36.7°.
To discuss the constant evaporation during stage 1, we show the stage 1 mass loss data for experiment 3 in Figure 12 in a higher resolution compared to Figure 11a. The experimental mass loss data are shown by black crosses. These data are smoothed by a Moving Average-filter with block size 10 using Mathematica 11. The smoothed data are shown by black circles. Linear regression of the smoothed experimental data yields a regression coefficient near 1 (= 0.99; black dashed line). Hence, a nearly perfect constant evaporation rate is shown in Figure 11b. Linear regression of the IP data (red points) yields also a regression coefficient near 1, that is, strict linear behavior of the IP-mass loss data and a constant evaporation rate, too.
As discussed by Mosthaf et al. (2014; see Figure 1 therein) one would expect a decreasing stage 1 evaporation rate for thin boundary layers (2 mm for experiment 3; see Table 2). This argument is based on a physical picture that lateral diffusion and interference effects between the diffusion sources (evaporating pores) are not Figure 11. (a) Mass loss of one effective square capillary versus time for experiment 3 (mean temperature 61°C, mean relative humidity (RH) = 40). The thin black solid line shows the best fit to the experimental data (blue dots). The red points represent the IP-mass loss data points. The blue straight line represents the best fit to the stage 1 data, and its slope determines the constant stage 1 evaporation rate. The theoretical mass loss curve (red curve, series expansion to third order; equation (11)) is compared to the analytical solution of the horizontal case (black curve, equations (14a) and (17a)). The inset shows the temperature and the relative humidity during the experiment. (b) Evaporation rate density versus time (log scale). The black curve is based on the best fit mass loss data and corresponds to the thin black line in (a). The red curve shows the sharp transition approximation of the single square capillary model. strong enough to establish a constant concentration profile along the surface. However, to quantify the interference effect, one has to consider the diffusion parameter, that is, the ratio BL thickness/distance of evaporating pores (see equation (3a) in Geistlinger & Leuther, 2018). If the diffusion parameter becomes larger than 1, then the lateral distance of diffusion sources (pores) has the same order of magnitude as the boundary layer thickness. At this threshold interference effects become dominant. For experiment 3 we have a diffusion parameter of 3.3. Hence, interference effects are strong enough to cause a constant steady-state concentration profile along the surface.

Detachment time t c
The experimental detachment time is 468 s ( Table 2). The theoretical t c value is given by equation (13), calculating the evaporation rate via equations (6a) and (6b) and the mass loss via equation (14a) for x 2 = 0. Solving the nonlinear equation (13) for the contact angle yields a contact angle of 36.5°. We emphasize that the x p and t c measurements are independent measurements. The equivalence of both values (deviation = 0.2°, relative error < 1%) demonstrate the consistency of the experimental measurements and, moreover, the validity of the θ-dependent SSC model. Both experimental values of 37°are strong hints that the SiO 2 surface causes partial wetting. Evaluating the evaporation experiments at lower temperatures in Exp. 1 (28°C) and Exp. 2 (42°C) yields similar contact angles: 37°for Exp. 2 and 39°for Exp. 1 and Exp. 4. The mean value over all contact angles listed in Table 1 is 37.9 ± 1.4°.
It is instructive to compare the theoretical contact angle derived from inverse modeling with the experimental contact angles of the flat surfaces of the Si-SiO 2 -wafer (48 ± 2°) and of the Borofloat-cover glass (23 ± 2°). Since the effective contact angle of our micromodel is the average of both contact angles, a simple arithmetic averaging yields: 3 × Si-SiO 2 -corners + 1 × glass corner)/4 = 42°. This yields a deviation of 10% between the effective experimental and the theoretical contact angle. This reasonable agreement supports the SSC model as an appropriate model for the description of the evaporation process in an ideal (no round corners!) square capillary during stage 1 evaporation.
We now discuss the temperature dependence of the stage 1 characteristics in order to test the consistency of the SSC model against the temperature variation. Therefore, we conducted evaporation experiments at three different temperatures: 28, 42, and 61°C. As a novel aspect, the presented theory takes into account both contact angle and temperature dependence of the water flux.
We first discuss the temperature dependence of the evaporation flux. Both the diffusion coefficient and the saturated vapor pressure increase with temperature and lead to a 2.5 larger evaporation flux density increasing the temperature from 28 (Exp. 1) to 42°C (Exp. 2). This is expected to cause a smaller width x p (T), as the width is inversely proportional to the evaporation rate (equations (12a) and (7c)). However, the x p values of the vertical experiments 1-3 are nearly constant and approximately 4 mm. Hence, the higher evaporation flux is compensated by a higher water flux. Obviously, the water flux will increase with temperature, because of the decreasing kinematic viscosity. However, this effect is not sufficient and increases the water flux only by a factor of 1.3. The main increase is caused by a decreased flow resistance factor α(θ) with decreasing contact angle from 39°to 37°. In this range α(θ) is a rather sensitive function of the contact angle (see Figure 6). The temperature-dependent decrease of the contact angle could be attributed to a decreased surface tension σ gw (T) with increasing temperature. Young's law, cos(θ) = σ w,s /σ gw (T), yields a smaller contact angle assuming that the strong water-solid interaction is temperature independent and neglecting the weak gas-solid interaction (Grant & Bachmann, 2002).

Time dependence of the mass loss curve
We now study the time dependence of the mass loss curve. We note that the experimental data shown in Figure 11 are renormalized data, that is, the mass loss and the volumetric evaporation rate for one effective capillary. The volumetric evaporation rate is given by the mass evaporation rate density j diff (Table 2) divided by N cap ρ w φ (φ-porosity, N cap -number of capillaries). There are two constraints to renormalize the original mass loss data: (i) the active surface area constraint and (ii) by the pore volume constraint. Setting the number of open surface pores to the number of capillaries (= 79), both methods lead to an effective capillary width. Note that the surface porosity (at the inlet) is not identical to the volume porosity of 0.46. The following discussion is based on the active surface area constraint, which gives an effective width of the square capillary of 0.3 mm.
The red curve in Figure 11a represents the theoretical mass loss curve for the vertical experiment 3 (Bo > 0). The mass loss equation (14a) is derived by using a series expansion third order for ξ(x) (equation (11)). Compared to the horizontal mass loss curve (thick black curve), the red curve fits the experimental data better. However, it still shows a significant deviation for longer times (>4,000 s).
This deviation indicates that during the later stage 2 period corner flow from small capillaries (pores) becomes important. The decreasing corner flow leads to a decreasing evaporation flux compared to the flux of the mean capillary of the SSC model. This presumably causes the SSC model to overestimate the evaporation flux. We tested this hypothesis by conducting a drainage experiment under the same capillary number. The resulting patterns were almost identical implying that the width of the PSD becomes increasingly important as evaporation proceeds.
We next analyze whether gravitational forces have an impact on the evaporation process, that is, causes a reduced mass loss. We therefore conducted a horizontal evaporation experiment under room temperature (Exp. 4) and compared it with the corresponding vertical experiment (Exp. 1; see Figure 13). The experimental mass loss under gravitation shows a significant reduced mass loss when compared to the horizontal experiment. The theoretical curve for the horizontal case (equation (14a)) shows reasonable agreement for t < 12,000 s. For later times the theoretical curve shows a significant deviation.
Based on the best fit of the mass loss data (thin black line in Figure 11a), we derive the evaporation flux density = (also denoted as evaporation rate) in millimeter per day shown by the thin black line in Figure 11b). Note that the experimental data show a smooth transition between stage 1 and stage 2 evaporation. The red curve in Figure 11b represents the sharp transition approximation of the SSC model.
An interesting question is whether the evaporation experiments exhibit a ffiffi t p behavior for the later stage 2 time interval. Figure 14 shows the mass loss curves versus ffiffi t p . All experiments are in excellent agreement with ffiffi t p behavior (regression coefficient R 2 = 1.0). This square root time dependence is usually obtained by infiltration theory (Brutsaert & Chen, 1995;Parlange et al., 1985), that is, by solving the nonlinear partial differential equation for the volumetric water content θ w (x,t) (see equations (26)-(28) in Geistlinger & Leuther, 2018). The solution gives with the desorptivity S d as a functional of the hydraulic functions of the porous media. The water flux at the surface decreases, because the water front recedes into the porous media, and hence, the flow has to overcome viscous forces over a longer flow distance. Hence, desorption theory considers (i) only the water movement, that is, diffusive transport within the gas phase is not considered, and (ii) there is no sharp evaporation front receding into the porous media. That means that the theory gives a continuous distribution of the volumetric water content between x 1 (percolation front) and x = 0 (surface).
The SSC model couples the water flux (water phase) with the vapor flux (diffusion within the gas phase), that is, it takes into account the receding water front x i (t) and the x i (t) dependence of the concentration gradient. It is interesting that the steady-state approximation of the SSC model yields

10.1029/2018WR024647
Water Resources Research the ffiffi t p dependence. We previously observed this ffiffi t p dependence for two real soils with high accuracy (relative error <2%; Geistlinger & Leuther, 2018).
One could argue that the steady-state mass flux (maximal mass flux for a given concentration gradient) can explain why the theoretical mass loss is too large compared to the experimental one. However, the velocity of the receding percolation front (Exp. 3: 2.5 × 10 −3 mm/s) is by 1 order of magnitude slower than the diffusion velocity of the water molecules in the gas phase (Exp. 3: 2.6 × 10 −2 mm/s). Hence, the diffusion process always adjusts the vapor concentration in the corner flow region to the new steady state. Therefore, it seems that the deviation is a consequence of the "stiff" solution of the effective capillary model (see discussion below).

Geometric characteristics of the stage 2 evaporation process
It is clear that geometric heterogeneity of the percolation front and the evaporation front will lead to heterogeneous diffusion lengths across the corner flow region. As the fast diffusion process within the gaseous phase will even this pore scale heterogeneity through lateral diffusion, an effective capillary model should be able to describe a diffusion-averaged quantity like the mass loss. The interesting question is whether the SSC model can describe averaged geometric characteristics of the evaporation process, like the mean evaporation and percolation front. Figure 15 compares the theoretical curves (red curve) with the experimental data (blue dots) for the stage 2 time interval.
Both the x i and the x p curves show strong deviations (up to 100% at t = 12,000 s), whereas the x i curve gives too small values and the x p curve too large values. Note that there is no free adjustable parameter, because the contact angle of 38°is determined by stage 1 evaporation characteristics.
The reason for these large deviations of the x i and x p values is the "stiffness" of the SSC solution illustrated in Figure 7. The model needs a too large flow length, that is, a too large width of the corner flow region, because it describes a tortuous corner flow path (shown in Figure 10) by a straight or stiff flow path of an effective capillary. However, the path has to be longer in order to be able to account for the same water pressure gradient. To understand why the model yields too small x i (t) values and too large x p (t) values, one has to consider the flux condition at the evaporation front (equation (9)). A decreasing evaporation rate (compare red and blue curves in Figure 7b) causes a decreasing slope of the solution ϕ(x) with time (compare red and blue curves in Figure 7a). This decreasing slope of the linear ϕ(x) results in a large change of the percolation front Δx 1 for a small change of the evaporation front Δx 2 . Therefore, the ratio Δx 1 /Δx 2 is always larger than 1 and increases with time. Hence, a small change of the evaporation front leads to a large increase of the mass loss (see Figure 7a). This forces the model to too small x i (t) values and to too large x p (t) values during stage 2. A realistic tortuous ϕ(x) solution is shown by the black dashed line in the ϕ(x) diagram of Figure 7a. This hypothetical ϕ(x) solution represents a tortuous path of the corner flow through a regular square lattice (see Figure 10). The slope satisfies the mass flux boundary condition (equation (9)) at x 2 (same slope as the linear solution) and leads to the same Δm w increase but to a larger x 2 value.

Summary and Conclusion
We conducted a series of evaporation experiments at different temperatures (23 to 61°C) in a micromodel (regular square 80 × 80 lattice with lognormal-distributed throats), which was produced by a new interval-based ICP-DRIE technology. This anisotropic etching technology guarantees high edge steepness (sharp corners) and a true mapping of the lattice structure with depth.
We visualized the tortuous corner flow paths between the percolation and evaporation front and analyzed the fluid-fluid patterns during the evaporation process. The mass loss was quantified both gravimetrically and by image processing.
Based on the steady-state assumption, we developed a θand T-dependent analytical theory considering corner flow within an effective square capillary and coupled it to the evaporation flux. Based on this theory, we Figure 15. Comparison between the power series solution (equation (11); red curves) and the experimental data (blue dots) for the mean evaporation front x i and the mean percolation front x p . The black curves represent the best fits of the analytical solutions (equations (17a) and (18)) to the experimental data.
were able to study the combined effect of contact angle and temperature on the evaporation process. We derived a novel explicit analytical power series solution of the Laplace equation that describes the steady-state corner flow in arbitrary direction. For the horizontal case (Bo = 0) described by the zero-order term of the solution, we obtained analytical expressions for the time dependence (i) of the mass loss (linear dependence for stage 1 and square root dependence for stage 2), (ii) of the mean evaporation front, and (iii) of the mean percolation front. For the vertical case and for Bo ≤ 10 −3 the power series solution shows fast convergence. Since our micromodel exhibits a relatively large effective pore diameter of about 0.4 mm, a third-order power series solution should be valid for most of the natural porous media, like soils, sands, sandstones, and fractured rocks.
We note that the general SSC model developed by Yiotis et al. (2012a) presents the results in terms of the hypergeometric function that has to be evaluated numerically to obtain the time-dependent characteristics of the evaporation process.
Our main conclusions are as follows: 1. Contact angle and temperature dependence The SSC model was able to describe the mass loss time behavior during stage 1 and stage 2 at different temperatures between 23 and 61°C. Interesting is that the model gives highly consistent contact angles derived from stage 1 characteristics: extension of the spanning water cluster and the critical detachment time. The experimental value of 38°± 1°is a strong indicator that the SiO 2 surface causes partial wetting. The consistent contact angles derived from stage 1 characteristics (deviation < 1%) and the reasonable agreement between effective experimental and theoretical contact angle (deviation = 10%) support the SSC model as an appropriate model describing the evaporation process in an ideal (no round corners) square capillary for stage 1 evaporation.
Any model that does not take into account the contact angle dependence will fail to explain the experimental x p and t c values or has to include another fitting parameter into the model, like the roundness factor.
For the early stage 2 time interval, the SSC model gives reasonable agreement between experimental and theoretical mass loss curves both for horizontal and for vertical evaporation. For the later-tage 2 time interval there is a significant deviation. An interesting experimental result is that all experimental mass loss curves show a ffiffi t p behavior in the later stage 2 time interval with high reliability (R 2 = 1.0). This ffiffi t p behavior is described well by the analytical solution, if used as a fit function.

Impact of gravitational forces: Horizontal versus vertical evaporation
To analyze whether gravitational forces have an impact on the evaporation process, that is, cause a reduced mass loss, we conducted a horizontal evaporation experiment under room temperature and compared it with the corresponding vertical experiment. The experimental mass loss under gravitation showed a significant reduced mass loss compared to that of the horizontal experiment. To the best of our knowledge, this is the first experimental demonstration that gravitation leads to a reduced mass loss during evaporation in 2-D micromodels.

Geometric characteristics of the evaporation process
The SSC model fails to describe the geometric characteristics of the evaporation process, like the mean evaporation and percolation front. Both the x i and the x p values show strong deviations during stage 2 (up to 100%), because the model describes a tortuous corner flow path by a straight flow path of an effective capillary.

Semiphenomenological models: IPE model versus SSC model
Our experimental study demonstrated that the SSC model gives good agreement between experimental and theoretical results during stage 1 and early stage 2 behavior. The SSC model is hereby based on renormalized mass loss data, that is, the effective width of the square capillary is determined by the active surface constraint. Consequently, observed agreement confirms the fundamental hypotheses of the SSC model. Namely, that the steady-state corner flow within a square capillary of effective width (= mean pore size and not the width of the PSD; hypotheses of the IPE model) controls the main characteristics of the evaporation process during stage 1 and early stage 2 behavior. These characteristics include extension of the spanning water cluster, detachment time, constant evaporation rate during stage 1, and 1/ ffiffi t p falling evaporation rate.
For the later stage 2 period the SSC model overestimates the evaporation flux, because evaporation from small capillaries becomes important implying that the hydraulic connected region of the PSD controls this evaporation period. The width of the effective capillary overestimates the size of small capillaries and hence the corner flow and the evaporation flux, too.
It is difficult to recommend which semiphenomenological model should be used to describe evaporation in real soils. One could use a certain section of the experimental PSD (= characteristic width of the PSD = hydraulic connected region) shown in Figure 4. The IPE model can be used to derive the stage 1 evaporation rate. We have shown before that there is some nonuniqueness (arbitrariness) choosing this characteristic PSD section; that is, different characteristic sections will lead to the same theoretical evaporation rate (Geistlinger & Leuther, 2018). The main drawback of the IPE model is that it fails to describe stage 2 and any geometric characteristics of the evaporation process. Based on these arguments, we would recommend the use of the SSC model. However, our evaporation study was based on a regular square lattice with ideal rectangular channel geometry. The mapping of this porous medium to an effective square capillary is straightforward, and we expect good agreement between theory and experiment. However, it is an unresolved problem how to map real porous media (like glass beads, sands, and soils) to an effective square capillary.
An alternative is a pore network model that considers the drainage of water within the pore channel and the buildup of water-filled corners and continuous tortuous corner flow paths. However, it seems to be unrealistic to develop such a pore network model that takes into account corner flow instability (snap-off) of about 5,000 grain corners within the corner flow region at each time step, to classify the corner flow into continuous and discontinuous corner flow paths ( Figure 10, dashed line section). Nevertheless, it would be an interesting question, if such an extended pore network model could describe the geometric characteristics of our experimental evaporation experiments.