A Dynamic and Fully Implicit Non‐Isothermal, Two‐Phase, Two‐Component Pore‐Network Model Coupled to Single‐Phase Free Flow for the Pore‐Scale Description of Evaporation Processes

We couple a dynamic, fully implicit non‐isothermal two‐phase, two‐component pore‐network model (PNM) to a free‐flow domain where the single‐phase flow Navier‐Stokes equations with component and energy transport are solved. A monolithic coupling scheme guarantees the conservation of mass, momentum and energy across the coupling interface and allows solving the system without any coupling iterations. We present a numerical example and show that the model is able to temporally resolve the pore‐local flow and transport dynamics of an evaporation process.


Introduction
Coupled system of free flow above a permeable medium appear in a wide range of technical and environmental applications. Prominent examples including two-phase flow range from fuel-cell water management (Gurau & Mann, 2009) and food drying (Defraeye et al., 2016;Verboven et al., 2006) to soil evaporation Or et al., 2013;Vanderborght et al., 2017) and salinization (Jambhekar et al., 2016;Mejri et al., 2017;Shokri-Kuehni et al., 2017). Modeling these types of systems is challenging due to the variety of length and time scales involved and the high physical complexity of the processes. Very often, the high computational demand of a discrete pore-scale description of the flow and transport phenomena within the porous domain limits numerical simulations to rather small domains. On the other hand, more efficient averaged approaches based on the concept of a representative elementary volume (REV) (Whitaker, 1999) may not provide sufficient pore-scale insight which might be required to fully resolve and understand the process. REV-type models may fail to account for sub-scale complexities such as pore-local water saturation distribution patterns which can have an influence on the total evaporation rates (Mosthaf et al., 2014;Shahraeeni et al., 2012). Pore-network models (PNM) (Fatt, 1956a(Fatt, , 1956b(Fatt, , 1956c capture pore-scale physics in an efficient yet relatively accurate manner (e.g. Oostrom et al., 2016) by considering a simplified but equivalent pore geometry. Using a quasi-static approach based on an invasion-percolation algorithm, they have been used to model the drying of porous media (see Prat (2011) for an excellent review). Here, viscous effects are neglected and the temporal dynamics of the process beyond capillary equilibrium remain unresolved. Based on the work of Thompson (2002), Joekar-Niasar and Hassanizadeh (2011) developed a dynamic two-phase PNM accounting for dynamics of the phase displacement process. Qin (2015) extended this concept in order to account for vapor transport. Beyhaghi et al. (2016) developed an iterative approach to couple a quasi-static PNM to a region of free flow. In previous work, we coupled a (Navier-) Stokes model to a single-phase PNM with component transport using a monolithic approach which does not require any coupling iterations  and extended the concept in order to account for pore-scale slip at the interface .
Here, we develop a dynamic, fully implicit, non-isothermal, two-phase, two-component PNM and couple it monolithically to a Navier-Stokes model. A numerical example of a PNM attached to a free-flow channel shows to temporal dynamics of the drying process where the pore-local influence of pressure-driven advective gas fluxes has a major impact on the overall system behavior.

Conceptual and Numerical Model
A two-phase, two-component pore-network model (describing the porous domain Ω PNM ) is coupled to a single-phase Navier-Stokes model (gas, used for the free-flow channel Ω FF ). Both models are non-isothermal. Water is the main component of the liquid wetting phase w, while the pseudocomponent air constitutes the main part of the non-wetting gaseous phase g. We consider local thermodynamic equilibrium (Helmig, 1997). The relevant constitutive relations and physical properties are implemented in the H2OAir fluidsystem of DuMu x Lauser, 2012), using data provided by the industry standard IAPWS (Wagner & Pruß, 2002) and Reid et al. (1987). This means that the phase densities, viscosities, enthalpies, etc., are solution dependent with respect to pressure, temperature and phase composition.
We use a fully implicit backward Euler scheme for the temporal discretization of the coupled model.

Momentum Balance
Assuming a compressible Newtonian fluid and neglecting the influence of gravity and dilation (Truckenbrodt, 1996), the free flow in sub-domain Ω FF is governed by the Navier-Stokes equations: Here, ϱ g and μ g = ϱ g ν g are the gas phase mass density and dynamic viscosity while v g and p g are the phase velocity and pressure. f represents an optional volume force.

Mole Balance
One molar balance equation is solved for each component κ ∈ {water, air}: x is the mole fraction of component κ within the gas phase and ϱ mol, g is the molar density of the latter. q κ is a molar sink or source term while  ,FF diff, mol, g j is the diffusive molar flux, which is approximated by Fick's first law here: Here, D g is the binary diffusion coefficient for water and air within the gas phase and  g X is the mass fraction of component κ. Note that the mass-based form of Fick's first law is used because Equation 1 yields a mass-averaged bulk velocity (Taylor & Krishna, 1993). Therefore, the diffusive mass flux needs to be converted into a molar flux by division by the molar mass M κ . We require   water air g g 1 x x .

Energy Balance
Non-isothermal systems require the formulation of an energy balance: The storage term includes the solution-dependent specific internal energy u g while the advective second term considers the specific phase enthalpy, takes into account the molar diffusive energy transport.
The third term describes the conductive energy flux via Fourier's law, where λ g is the phase heat conductivity and T the phase temperature. q e is an energy sink or source term.

Numerical Model
The free-flow model is discretized in space using a staggered-grid finite volume approach, also known as MAC scheme (Harlow & Welch, 1965), which provides inherently stable and oscillation-free solutions without the need of any stabilization techniques (Versteeg & Malalasekera, 2007). Scalar quantities like pressure, density or mole fractions are stored at the centers of the control volumes C of the primal grid , while the velocity degrees of freedom are located at the primary cells' faces around which dual-grid ( *  ) control volumes * { , , } x y z C are constructed. We refer to (Schneider et al., 2020) for details on the discretized form of Equation 1. The primary variables of the free-flow model are p g , v g , T, and water g x .

Pore-network Model for Ω PNM
The physics considered in Ω PNM correspond to the REV-scale model of Mosthaf et al. (2011), adapted to a pore-scale description of the system: • two phases α ∈ {wetting liquid (w), gas (g)}, two components κ ∈ {water, air} • miscible, non-isothermal system with phase change • creeping (Re < 1) flow, no influence of gravity The dynamic model is based on the work of Joekar-Niasar et al. (2010) and Qin (2015) with some important modifications and extensions. Instead of a semi-implicit approach, we use a fully implicit, fully coupled model formulation and consider individual balance equations for both components and energy.
We consider cubic pore bodies and throats with rectangular cross-section in this work, but the model can be easily adapted for other geometries, too.

Momentum Balance
The (fully developed) volume flow of phase α ∈{l, g} through pore throat ij connecting pore bodies i and j (Joekar-Niasar et al., 2010) is given by The phase pressures defined at the pore-body centers are given by p α,i and p α,j while g α,ij is a phase-specific pore-throat conductance factor which depends on the throat geometry, the fluid viscosities and the phase distribution (i.e., throat-local capillary pressure).
If only one phase is present, g α,ij (Patzek and Silin, 2001) is given by with k = 0.5623 and the shape factor G = 1/16. The throat's total cross-sectional area and its length are given by A tot,ij and l ij . The pressure drop within the pore body is neglected. If both phases are present, the wetting layer conductance is given by Ransohoff and Radke (1988), where r am = γ/p c is the radius of curvature of the arc menisci, A w,Λ the cross-sectional area of the wetting layer in corner Λ, and ξ Λ is a dimensionless resistance factor (Zhou et al., 1997). The gas-phase conductance at the center of the throat (Bakke and Øren, 1997) is A g,ij is the cross-sectional area of the gas phase and r ij is the throat's inscribed radius. Before the gas phase can enter a water-filled pore throat, the entry capillary pressure (Mayer & Stowe, 1965;Princen, 1969aPrincen, , 1969bPrincen, , 1970) needs to be overcome. We assume a contact angle θ = 0 and an interface tension of γ = 0.0725 kg s −2 .
Snap-off (Blunt, 2017) occurs if the capillary pressure falls below where β = π/4 is the corner half-angle.

Mole Balance
A component balance is formulated for each component κ: S α,i is the saturation of phase α at pore body i. The diffusive fluxes are again described by Fick's first law: The constitutive relations mentioned above close the system of equations, that is, phase transfer processes such as evaporation or condensation are not modeled explicitly (Nuske et al., 2014), but implicitly accounted for under the assumption of a local thermodynamic equilibrium per pore body (Helmig, 1997).
Numerical problems arise if one of the phases vanishes, for example, when the drying of soil is considered, as the two-phase system degenerates and S α and   x of the vanished phase lose their physical meaning. A possible solution strategy is switching the local set of primary variables such that they remain physically meaningful .
Such models can suffer from decreased numerical convergence behavior and oscillations if the switch is triggered and the number of present phases is about to change locally .
In this work, we follow a more stable but computationally more involved approach solving for an additional constraint for each phase α (Lauser, 2014;Lauser et al., 2011): This NCP (non-linear complementarity problem) model has a fixed set of primary variables which is p g , S w and the fugacities f κ of all components κ.
WEISHAUPT AND HELMIG

10.1029/2020WR028772
Equation 15 is solved in terms of in this work, but other conditions, such as the Fischer-Burmeister (Fischer, 1992) function, fulfill Equation 15, too.

Energy Balance
The energy balance is given by: The molecular diffusive energy transport is taken into account by  h . In analogy to Equation 5, the specific phase enthalpies depend on the composition: Note that only energy fluxes within the pore space (the fluid) are balanced here.

Closure Relations
We use the pore-local p c − S w relation for cubic pore bodies as given by Joekar-Niasar et al. (2010): Finally, we require S w,i + S g,i = 1 and     water air , ,

Numerical Model
We consider a dynamic and fully coupled approach such that both saturations and pressures are solved for simultaneously, using the Newton-Raphson method. The time step size is chosen adaptively based on the convergence behavior of the non-linear solver (e.g., Baber et al., 2012).
After each Newton iteration, the invasion state of each throat is checked. If the throat was not invaded in the previous time step and if one of the capillary pressures of the adjacent pore bodies is greater than the entry capillary pressure p c,e (see Equation 11), the throat is marked as invaded and the set of constitutive laws for describing fluid flow is adapted accordingly. If an invasion or snap-off event occurs, we force the Newton scheme to do at least one additional iteration before the solution is considered to be converged. In order to increase numerical robustness we regularize Equation 19 linearly for very high saturations and we modify the curve such that p c,i (S w,i = 1) = 0. If the change of saturation within a pore body between two consecutive time steps is larger than a given threshold (e.g., 30%), the time step will be rejected and repeated with a reduced time step size.
During preliminary numerical studies, the Newton algorithm failed to converge for S w approaching zero. To remedy this, we also introduced a primary variable switch for the NCP model, replacing p g by p w as primary variable if S w falls below a given threshold and vice versa. We found that this very simple mechanism does not impair the numerical robustness of our model, in contrast to above-mentioned primary variable switch models where initial saturations of the appearing phases need to be guessed rather arbitrarily in case of a switch event . For sake of simplicity, we will continue to use the term NCP model even though the modified version of it is meant.

Coupling Conditions
Appropriate coupling conditions ensure the continuity of mass, momentum and energy at the interface between the free flow and the pore-network model (Hassanizadeh & Gray, 1989).
All coupling conditions are formulated for creeping (Re < 1), laminar flow at the interface. We follow largely  with certain adaptions made to account for the pore-scale coupling. The superscripts FF and PNM refer to the interfacial quantities of the free-flow domain and the pore-network model respectively. n is a unit vector normal to the coupling interface, pointing outside the own model domain.
REV-scale coupled models (e.g. Layton et al., 2002;Mosthaf et al., 2011) consider the exchange of mass, momentum and energy across the entire coupling interface Γ, that is, the common surface between the free-flow and porous-medium model. In contrast to that, the free-flow and the pore-network model are only coupled at discrete locations, the pore-scale coupling interfaces Γ i which are associated with each pore body i intersecting with the boundary of the free-flow domain (Figure 1).
These pore bodies are assumed to be cut in half by the interface and only the interior part of the volume is considered. We assume that the coordination number of pore bodies connected to the free-flow domain is always one, that is, only one pore throat is connected to them. The free-flow grid cell sizes are chosen such that they are perfectly aligned with Γ i , that is, they are flush with the borders of the interface pore. Any number >0 of free-flow cells per boundary pore may be chosen. At the location of solid grains (no intersecting pore body), a no-slip/no-flow condition for the free flow is assumed, as indicated for the gray solid part in Figure 1b.

Momentum Balance
The momentum interface conditions are explained in detail in Weishaupt et al. (2020). We require mechanical equilibrium at the interface and assume creeping flow (Re < 1), yielding for the normal force component. The tangential momentum transfer from the pore-network model to the free-flow domain is described by an approximation similar to the widely used Beavers-Joseph interface slip condition for REV-scale models (Beavers & Joseph, 1967), adapted for a pore-local description of the slip .
WEISHAUPT AND HELMIG 10.1029/2020WR028772 6 of 14 is the basis of the interface's tangent plane. The tangential component of the pore-body interface velocity is approximated as Q g,ij is the gas phase volume flow through pore throat ij while |Γ i | is the area of the discrete coupling interface. n ij is a unit normal vector parallel to the throat's central axis and pointing toward the interface. Note that this is a simplification which does not take into account potential deflection effects of the fluid flow leaving the pore throat and entering the pore body. We furthermore assume that the gas phase completely occupies Γ i , neglecting potential two-phase effects directly at the interface. Note that Equation 23 does not impair mass conservation as it is only used for the approximation of the tangential momentum transfer. The material parameter β pore is determined numerically in a preprocessing step .

Chemical Equilibrium and Mole Flux Conservation
For compositional flow, we require the conservation of the individual component fluxes across the interface for each component κ: We assume local chemical equilibrium ) within a pore body and at the interface, which serves as a Dirichlet-type coupling condition for the free-flow model.
The diffusive molecular flux  diff, mol, g j is approximated using the mass fraction gradient within the free-flow domain directly at the interface: Here, the superscript FF + 1 denotes the quantities of the first free-flow grid cell adjacent to the interface, while 0.5[Δh] FF+1 is the distance between the interface and said grid cell's center. In the presence of multiple phases within the pore body at the interface we assume that all phase-transfer processes occur within the pore body under local thermodynamic equilibrium. This means that liquid water vapourizes in the pore body and only the diffusive and advective flux of water vapor within the gas phase is considered for the mass transfer between the pore-network model and the free-flow model. This limitation might become relevant in case of pronounced liquid-phase flow toward the interface driven by viscous forces (e.g., when an external liquid-phase pressure gradient is applied) which we do not consider in this work. The coupling condition could be extended in future work such that liquid water is allowed to leave the network. One could either employ a two-phase flow model concept for Ω FF or consider the formation of liquid drops (e.g. Ackermann et al., 2020;Baber et al., 2016) at the interface.

Energy Balance
For non-isothermal flow, we require the conservation of energy fluxes across the interface: As discussed earlier in the context of diffusive mass exchange between the sub-models, we assume that the gas phase present in the free-flow domain fully occupies the coupling interface Γ i and therefore, all diffusive, conductive and convective energy transfer only occurs through phase g. The conductive heat transfer is calculated in analogy to Equation 26 as we demand local thermal equilibrium at the coupling interface:

Implementation and Software
The model in implemented in DuMu x Heck et al., 2019;Koch et al., 2020), a free and open-source toolbox for porous-media applications based on Dune (Bastian, Blatt, Dedner, Engwer, Klöfkorn, Kornhuber, & Sander, 2008;Bastian, Blatt, Dedner, Engwer, Klöfkorn, Ohlberger, & Sander, 2008;Sander et al., 2017). We follow a monolithic coupling concept such that all balance equations are assembled into one linearized system of equations which is solved by UMFPack (multifrontal LU factorization, (Davis, 2004)). This means that no coupling iterations are required and the scheme inherently conserves mass, momentum and energy. Details on the system matrix can be found in Weishaupt et al. (2020).

Evaporation from a Network of Pores
We demonstrate the model's capabilities on the example of water evaporating from a porous domain Ω PNM into a free-flow channel Ω FF where a constant stream of dry and warm air is supplied from the channel inlet on the left (Figure 2a). The geometry is inspired by a single-phase flow microfluidic experiment (Terzis et al., 2019;Weishaupt et al., 2020). Further numerical examples may be found in Weishaupt (2020), were we also consider cases with a stagnant air atmosphere and throats of different shapes. The channel is 1.5 × 10 −2 m long, 0.5 × 10 −2 m wide and 4 × 10 −4 m high. The rather high aspect ratio of 12.5 between the channel's width and height allows the use of a quasi-3D approach where a two-dimensional numerical grid is employed to describe the three-dimensional physical domain, thus increasing the model's efficiency. Assuming a parabolic flow profile along the omitted height, Flekkøy et al. (1995) proposed a drag term (added to Equation 1) which accounts for the wall friction of the virtual frontal and rearward boundary: h Ω is the virtual height of the model domain while c is a constant which determines whether the maximum velocity at the central plane of the channel at 0.5 h Ω (c = 8) or the height-averaged one (c = 12) is recovered. We chose the latter here. This approach has been applied successfully for a number of different applications with thin, slit-like domains (Venturoli & Boek, 2006;Laleian et al., 2015;Kunz et al., 2015;Class et al., 2020). The channel is discretized with 30 cells in vertical and 133 cells in horizontal direction such that three grid cells coincide with one pore body at the interface. The pore-scale slip factor β pore = 19,134 for Equation 22 was determined numerically in a preprocessing step .
We apply fixed pressures on the left and the right side of the channel while the top and bottom of the domain (whit the exception of Γ i ) feature a no-flow/no-slip boundary condition for mass and energy (Figure 2a). An outflow condition for mass and energy holds at the right boundary.
We chose initial and boundary values of v g, init = 0, p g, init = p g, out = 1 × 10 5 Pa,   water water g, init g, in x x 4.7 × 10 −3 (corresponding to a relative humidity of 0.2) and T init = T in = 293.15 K. Air flow is induced by setting p g, in = p g, out + 1 Pa, yielding a maximum velocity of around 5.1 × 10 −2 m/s and Re ≈ 17 in Ω FF .
The solid matrix of the porous medium is assumed to be perfectly insulating and thus no heat exchange between the fluids in the network and the solid matrix (which is not modeled anyway) takes place. This process needs to be addressed in future work (e.g., Surasani et al., 2008).
All sides of Ω PNM are closed and no-flow conditions for mass and energy apply. Initially, Ω PNM is fully water-saturated (S w, init = 1,  air w, init x 1 × 10 −5 ), except for the pores at the coupling interface Γ FF where S w, init = 0.99 such that there is already some gaseous air in the pore bodies which is required conceptually for the coupling conditions described in Section 2.3. p g, init = 1 × 10 5 Pa and T init = 293.15 K are chosen as initial conditions for the pressure and the temperature. As explained above, the coupling conditions do not permit liquid water to leave the network. This would only be a conceptual limitation for fully water-saturated pores at the coupling interface which never occurred during the simulation run. Figure 2 shows the temporal development of the evaporation process and indicates the formation of a vapor boundary layer above the interface between Ω PNM and Ω FF . As seen in Figure 2c  and ultimately vanishes as the medium dries out completely, beginning from the upstream (left) region of the interface.
The invasion of the air phase into the network starts on the left side of the interface, where the air humidity in Ω FF is lowest and the largest pore-local evaporation rates occur. The air phase spreads vertically in Ω PNM until it reaches the network's bottom and continues to invade throats in a rightward direction. It temporally encloses a patch of completely saturated pores, as shown in Figure 2b, before the network eventually dries out. This is probably a boundary effect related to the decreased coordination number of the lower-bottom pores. If an interior pore starts getting drained from one throat, there are three possible neighbor pores which could replenish the exiting mass of water. For boundary pores, only two such potential neighbor pores exists such that the chance for effectively decreasing the wetting-phase saturation and thus increasing p c above p c,e is higher.
The global evaporation rate (i.e., the sum of all pore-local evaporative fluxes) is shown by the red solid curve in Figure 2d. For comparison, we also plotted the corresponding data for the same setup with a stagnant air atmosphere (closed inlet and outlet, open top with  water g, top x 4.7 × 10 −3 , T top = 293.15 K and p g, top = 1 × 10 5 Pa) as discussed in detail in Weishaupt (2020). The continuous supply of warm, dry of air in Ω FF leads to a higher evaporation rate for the case with air flow, resulting in only around 14 h needed to completely dry the medium compared to 25 h for the resting atmosphere case. The latter features a rather typical evaporation rate over time: the initial fall of the curve due to temperature equilibration at the interface is followed by a long period of constant stage-I evaporation (Or et al., 2013;Yiotis et al., 2007) due to the presence of liquid water at the surface which only vanishes at the very end of the simulation. In contrast to that, the temporal development of the evaporation rate for flowing air may seem unexpected at first glance: after the initial fall and a short constant phase, the evaporation rate actually increases again over time. This behavior results from chosen geometry and boundary conditions, inspired by a microfluidic experiment. We consider a rather coarse (and hence highly permeable) porous medium which is furthermore constrained by walls on all sides but the coupling interface. As shown experimentally (Terzis et al., 2019) and numerically  for single-phase flow, this results in a distinct pressure-driven U-shaped flow pattern in the porous domain, where the fluid vertically enters Ω PNM on the left side and leaves it on the right side. This distinctive pattern can also be observed here once the gas phase is present in Ω PNM . Figure 3 shows the local advective flow of the gas phase Q g,ij within the network's throats (A-C) as well as the pore-local evaporation rates at the interface for different times (D). Initially, only the vertical throats at the interface are invaded at t = 1 h (no advective gas flux in Ω PNM , not shown) which corresponds to the expected distribution of local evaporation rates along the interface shown by the blue curve in Figure 3d: Evaporation is strongest at the leftmost pore, where the air in Ω FF is driest, and continues to decrease along the interface as the air gets more and more saturated with water vapor. The situation changes, however, as soon as the first row of horizontal throats next to the interface gets invaded (see Figure 3a, t = 2 h) which establishes a flow path for the air phase. The arrows indicate that air now enters the network at the leftmost throat #1, flows parallel to the interface and exits the network again at the rightmost throat #19. This strongly increases the evaporative flux at throat #19, shown by the orange curve in Figure 3d. The flow of air pushes the water vapor out of Ω PNM , and hence the evaporation process is no longer only diffusion-driven but dominated by advection. This becomes even more clear as the air phase continuous to invade more and more throats, creating new pathways for advective gas flow in the network, as shown by Figure 3b (t = 9 h) and the green line in Figure 3d.
The latter exhibits a drop at pores #14-#17 due to the locally decreased outflow of gas which is caused by the given saturation distribution in the network. When all throats are invaded by the gas phase (Figure 3c, t = 11 h), the flow field in Ω PNM essentially resembles the one described for single-phase flow in Weishaupt et al. (2019): the air enters the network along the left half of the interface and leaves Ω PNM symmetrically across the right half of the interface, which explains the strongly elevated vapor fluxes on this side, as shown by the red curve in Figure 3d. Diffusion-driven evaporation still takes place due to the presence of liquid water at the interface. At the left half of the interface, diffusion and advection are actually two competing mechanisms as the vapor concentration gradient yields a diffusive flux toward Ω FF , reversely to the inbound advective flow. The red and green curves in Figure 3d also show this reduction in evaporation rates for pores #1-#10. These pore-scale effects could not have been observed using standard quasi-static pore-network models (Beyhaghi et al., 2016) which do not resolve fluid flow dynamics beyond capillary equilibrium states.
In summary, we could observe the formation of a water vapor boundary layer above the interface. Evaporation happens significantly faster compared to a setup with a resting atmosphere. The invasion of the gas phase into the network occurs non-uniformly, creating pathways for advective gas flux in Ω PNM . This leads to a rather atypical behavior of the evaporation rate over time. The presence of liquid water at the surface technically qualifies the process as stage-I evaporation where constant or falling evaporation rates would be expected (Or et al., 2013). However, the evaporation rate actually increases over time due to advective water vapor transport out of Ω PNM which is a result of the chosen geometry and boundary conditions. Evaporation is no longer only diffusion-controlled but governed by advection, too.
While not shown here, the novel model was also applied to a case with stagnant air and throats with round cross-sectional shapes. This considerably decreased the global evaporation rate due to the absence of throat corner flow (which is only possible for angular throats) such that mass transfer was soon limited by diffusion which corresponds to stage-II evaporation (Or et al., 2013). We furthermore compared the NCP approach presented here (Equation 15) with a primary variable switch method  and found that the latter either fails to converge at all for a complete dry-out of the system for certain configurations or needs more than twice the CPU time required by the NCP model. Finally, we investigated a case with stagnant air and a heterogeneous pore network. Here we could observe a strong decrease in time step sizes for throat invasion events which limited the model's overall performance. We will address this in future work. We refer to Weishaupt (2020)

Conclusions
We have presented a coupled model concept for simulating evaporation processes on the pore-scale. A freeflow domain, governed by the Navier-Stokes equations, is coupled to a dynamic, fully implicit two-phase, two-component pore-network model. A numerical example illustrates the model's capability to resolve the temporal dynamics of the drying of a network. As given by the setup and the boundary conditions, the process strongly influenced by advection, once the air phase was able to invade the network and create continuous pathways. This application showed to importance of considering dynamic effects beyond capillary equilibrium.
The model can be a valuable tool to deepen the understanding of complex pore-scale physics in a rather efficient manner. It is especially suitable for planning and evaluating microfluidic experiments (mm to cm scale) and it could be an interesting option for fuel cell modeling. The pore-scale results and insights could be used to derive upscaled relations for REV-scale concepts allowing the simulation of larger domains. As a next step, we aim at improving the model's numerical performance for heterogeneous media, and we will use microfluidic evaporation experiments for validation purposes.

Data Availability Statement
The source code and instructions for reproducing the data presented in the last section is publicly available as dumux-pub module (https://git.iws.uni-stuttgart.de/dumux-pub/weishaupt2020a).