Drying of Nanoporous Media by Concurrent Drainage and Evaporation: A Pore Network Modeling Study Haiyi Wu,1 Chao Fang,1 Rui Wu,2 and Rui Qiao1,* 1 Department of Mechanical Engineering, Virginia Tech, Blacksburg, VA 24061, USA 2 School of Mechanical Engineering, Shanghai Jiao Tong University, Shanghai, 200240, China Abstract. Drainage and evaporation can occur simultaneously during the drying of porous media, but the interactions between these processes and their effects on drying are rarely studied. In this work, we develop a pore network model that considers drainage, evaporation, and rarefied multi-component gas transport in nanopores. Using this model, we investigate the drying of a liquid solvent-saturated nanoporous medium enabled by the flow of purge gas through it. Simulations show that drying progresses in three stages, and the solvent removal by drainage effects (evaporation effects) becomes increasingly weak (strong) as drying progresses through these stages. Interestingly, drainage can contribute considerably to solvent removal even after evaporation effects become very strong, especially when the applied pressure difference across the porous medium is low. We show that these phenomena are the results of the coupling between the drainage and evaporation effects and this coupling depends on the operating conditions and the stage of drying. Keywords: porous media; drying; pore network models; multiphase flows; multicomponent gas transport * To whom correspondence should be addressed. Email: ruiqiao@vt.edu. 1 I. Introduction Drying of porous media including building materials and pharmaceuticals is important in engineering applications ranging from material manufacturing and oil extraction to soil remediation. 1-3, 4 Although many technologies including thermal drying and infrared radiation drying have been developed,5 they often suffer from limitations such as high energy cost, low throughput, etc. Hence, there is a long-standing need to improve existing methods or to develop new technologies to overcome these limitations. Addressing this need through trial-and-error experiments is often costly and ineffective, and numerical modeling can be helpful. Modeling of the drying of porous media, however, is challenging because coupled multiphase heat and mass transfer at scales from nanometer to centimeters must be simulated.6-7 To tackle the multi-faceted (multiphase, multiscale, and multiphysics) drying processes in porous media, many numerical models have been developed. These models can be loosely classified into two categories: macrohomogeneous models and pore scale models. In macrohomogeneous models, the porous media are treated as a continuum with volume-averaged or homogenized properties. These models can deal with large scale problems easily. 8-9 However, because the microstructures of porous media and the physical processes in them (e.g., viscous fingering and corner flows) are not resolved but modeled heuristically using empirical or experimental sub-models, they offer limited insight into the fundamental physics of the drying process and often lack predictive power in novel situations. In pore scale models, the heat and mass transfer in the microscale geometry of porous media are considered. In some models, the original heat and mass transfer equations such as the Navier-Stokes equation are solved in porous media.10,11,12 In other models,13 porous media are represented using various building blocks and the transport processes are modeled based on rules derived from fundamental heat and mass transfer laws. Amount these models, pore network models have received much attention because they offer a good balance between capturing pore scale physics (hence drying) and being computationally efficient. In pore network models, the pore space is simplified as discrete elements consisting of pore bodies interconnected with each other by pore throats. The first pore network model for studying drying was built in 195413 and many improvements have been made since then. 14-19 Key to these models is the consideration of fluid and mass transfer at the pore scale and in the interconnected pore systems. For example, the transport of vapor obeys convection-diffusion equations and the continuity equation is applied to the liquid phase and non-condensable gas phase;14, 17 the pattern of receding liquid-vapor interfaces due to evaporation has been described using the invasion percolation (IP) model. 20-21 While existing pore network models have been successful in analyzing many porous media drying problems, new problems with unique features not explored previously continue to appear and the existing models need to be extended to study these problems. 2 In this work, we extend pore network models to study the drying of nanoporous media assisted by purge gas flowing through them. Figure 1 show a schematic of the drying problem studied here, which is part of the solvent recovery step in the dewatering-by-displacement technology.22 A bed of particles with diameter often smaller than 2 m is initially saturated with a volatile solvent (e.g., pentane) and a purge gas is driven through the particle bed to remove the solvent. At the beginning, the injected purge gas penetrates into the particle bed to displace the liquid solvent, just as in the conventional drainage process. Eventually, the purge gas breaks through the particle bed and gas transport pathways are formed across the porous matrix. After this, evaporation occurs within the porous matrix and the vaporized solvent is “flushed” downstream by the purge gas. Meanwhile, the purge gas continues to drive the remaining liquids out of the porous matrix. Although the extensive pore network modeling in the literature has led to useful fundamental understanding of porous media drying, they offer limited insights for this drying problem because of its several unique features: 1. Drying is assisted by the flow of purge gas through the porous media, which differs from the widely studied scenario where drying is assisted by gas streams flowing over the porous media’s surface; 2. The drainage of liquid solvents and evaporation of solvents can occur concurrently in the porous media and become tightly coupled; 3. Because of the small particle size, the pores between them can have diameter of a few hundred nanometers.23 Therefore, the gas slippage and Knudsen effects, rarely considered in previous drying research, can no longer be neglected in the gas transport process; 4. The liquid solvent is highly volatile and their vapor pressure can be comparable to that of the purge gas. The multicomponent gas transport behavior can involve strong couplings between the purge gas and vaporized solvents, which cannot be described well using the classical convection-diffusion equations. Figure 1. Drying of a nanoporous medium with purge gas flowing through it. (a) The porous medium, a particle bed, is initially saturated with a volatile liquid solvent. (b) The purge gas displaces liquid solvents from the porous medium, and the solvent evaporation is negligible. (c) The purge gas breaks through the porous medium and evaporation of the liquid solvents becomes important; drainage and evaporation can occur concurrently and become strongly coupled. 3 In this work, we study the drying of nanoporous media assisted by purge gas flowing through them using a pore network model. The rest of the manuscript is structured as follows: In Section II, we detail our model that accounts for the concurrent the drainage and evaporation processes during drying and the other unique features mentioned above. In Section III, we present the simulation results for drying, reporting the evolution of the drying rate and discussing the couplings between the drainage and evaporation effects. Finally, conclusions are presented in Section IV. II. Model Formulation To study the drying problem shown in Fig. 1, we assume that the nanoporous media is initially (𝑡 = 0) fully saturated with a volatile liquid solvent. The solvent liquids have a zero contact angle on solid surfaces and the purge gas has negligible solubility in the solvent liquids. Below, we first present the setup of the pore network, then discuss the physical and numerical models, and finally present the numerical algorithm. A. Pore network Since we focus on the coupling between drainage and evaporation, which is the most unique aspect of the dyring problem shown in Fig. 1, the geometry of the porous media is simplified as a two-dimensional (2D) rectangular block represented by a pore network. 24-25 The pore network lattice consists of pore bodies connected to their four neighbors via pore throats (see Fig. 2a). The upstream (downstream) boundary of the pore network is connected with a upstream (downstream) reservoir. Periodic boundary conditions are applied on the left and right boundaries. The pore bodies have cubic shape and the pore throats have circular cross-sections. A schematic of two pore bodies connected via a pore throat is presented in Fig. 2b. Following previous work,24-25 the size distribution of pore bodies and pore throats both follow a truncated log-normal distribution: 𝑓 (𝑟𝑖 ) = 𝑟 ln 𝑖 𝑟 √2 exp(−0.5( 𝜎𝑚 )) 𝑟 𝑟 √𝜋𝜎 2∙𝑟𝑖 (erf(ln 𝑚𝑎𝑥 ⁄√2𝜎 2)−erf(ln 𝑚𝑖𝑛 ⁄√2𝜎 2)) 𝑟𝑚 (1) 𝑟𝑚 where 𝑟𝑖 is the radius of the inscribed sphere of a cubic pore body i (see Fig. 2b); in other words, the length of the pore body is 2𝑟𝑖 . 𝜎 is the standard deviation. 𝑟𝑚 , 𝑟𝑚𝑖𝑛 and 𝑟𝑚𝑎𝑥 are the mean of the inscribed sphere radius, the lower bound and upper bound of the truncation, respectively. The spacing between adjacent layers in the pore network in the x- and y-directions are denoted as ∆𝑆𝑥,𝑖 and ∆𝑆𝑦,𝑖 , respectively: ∆𝑆𝑥,𝑖 = 𝑎𝑥 max{𝑟(𝑖, 𝑗) + 𝑟(𝑖 + 1, 𝑗), for 𝑗 = 1: 𝑛𝑦 } , 𝑖 = 1: 𝑛𝑥 ∆𝑆𝑦,𝑗 = 𝑎𝑦 max{𝑟(𝑖, 𝑗) + 𝑟(𝑖, 𝑗 + 1), for 𝑖 = 1: 𝑛𝑥 } , 𝑗 = 1: 𝑛𝑦 (2𝑎) (2𝑏) where 𝛼𝑥 and 𝛼𝑦 are prefactors set to 𝑎𝑥 = 𝑎𝑦 = 1.2 and 𝑛𝑥 (𝑛𝑦 ) is the number of pore bodies in the x- (y-) direction. The length of the pore throats is then determined via 4 { 𝑙𝑥,𝑖 = ∆𝑆𝑥,𝑖 − (𝑟(𝑖, 𝑗) + 𝑟(𝑖 + 1, 𝑗)) for 𝑗 = 1: 𝑛𝑦 , 𝑙𝑦,𝑗 = ∆𝑆𝑥,𝑗 − (𝑟(𝑖, 𝑗) + 𝑟(𝑖, 𝑗 + 1)) for 𝑖 = 1: 𝑛𝑥 , 𝑖 = 1: 𝑛𝑥 𝑗 = 1: 𝑛𝑦 (3𝑎) (3𝑏) Figure 2. The pore network model. (a) A schematic of the pore network. (b) A sketch of two pore bodies connected by a pore throat. (c) The pore body size distribution. (d) The pore throat size distribution. B. Physical and numerical models Below, we present the models for the drying of nanoporous media assisted by purge gas flowing through them. Key assumptions for model development, the models for liquid and vapor transport, the models for liquid vaporization, initial and boundary conditions, rules for events in the pore network model, and the computational algorithms are discussed. B.1 Assumptions The following assumptions are made in our pore network model: 1. Pore bodies have finite volume but zero hydraulic resistance for liquid and gas transport. 2. Pore throats have negligible volume and their hydraulic resistance corresponds to fully developed flows in ducts with the same size. A pore throat has only two states: filled with liquid or filled with gas. Meniscus is not tracked inside pore throats and the time required for filling/emptying a single pore throat by liquid/gas is negligible. 3. Liquid solvents are incompressible and the solid phase is considered as a rigid body. 5 4. For the gas and liquid flows in pore throats, the inertia effects are negligible. 5. Gravity effects are negligible because the drying process occurs at microscale. 6. Drying occurs under isothermal conditions. 7. The vaporization of liquid solvents at the corner of a pore body can occur if more than two throats connected with that pore body is occupied by the gas phase. For pore bodies in which vaporization does not occur, the density of the solvent vapor is set to its saturation density at the same temperature. These assumptions have been adopted extensively in prior pore network simulations of porous media drying because they allow capture the essential physics of the drying process to be modeled with modest cost and computer memory.17, 21, 24, 26-27 B.2 Transport model We adopt the two-pressure algorithm24, 28-29 for solving the pressure field for both liquid and gas phases. Hereafter, variables corresponding to the liquid and gas phases are identified using superscripts “l” and “g”, respectively. Because the gas phase is generally a two-component mixture of the purge gas (species 1) and the solvent vapor (species 2), variables corresponding these components are identified using superscripts “g1” and “g2”, respectively. The local capillary pressures in the pore bodies are defined and approximated as24, 29 𝑔 𝑃𝑖𝑐 (𝑆𝑖𝑙 ) = 𝑃𝑖 − 𝑃𝑖𝑙 = 2𝛾 𝑟𝑖 (1−exp(−6.83𝑆𝑖𝑙 )) (4) where 𝑃𝑖𝑔 𝑃𝑖𝑙 , 𝑃𝑖𝑐 and 𝑆𝑖𝑙 are the gas phase pressure, liquid phase pressure, capillary pressure, and the saturation of the liquid phase in the pore body i. 𝛾 is the interfacial tension. Volume balance: A volumetric flux is assigned to a pore throat ij for both gas and liquid phases and a volume balance is required for each pore body i: 𝑔 𝑖 𝑙 ∑𝑁 (𝑄𝑖𝑗 + 𝑄𝑖𝑗 ) = 0 𝑗=1 (5) where 𝑁𝑖 is the number of pore bodies connected to the pore i, 𝑄𝑖𝑗𝑙 (𝑄𝑖𝑗𝑔 ) is the volumetric flux of liquid (gas) through the pore throats ij. Liquid transport: The volumetric flux of the liquid phase through a pore throat ij is 𝑙 𝑄𝑖𝑗 = −𝐾𝑖𝑗𝑙 (𝑃𝑖𝑙 − 𝑃𝑗𝑙 ) (6) where 𝑃𝑖𝑙 is the liquid pressure of pore body i and, 𝐾𝑖𝑗𝑙 is the conductivity of the pore throats filled with 𝜋𝑟4 liquid. 𝐾𝑖𝑗𝑙 = 8𝜇 𝑖𝑗𝑙 since we assume a Hagen-Poiseuille flow in pore throats. 𝜇𝑙 , 𝑙𝑖𝑗 , 𝑟𝑖𝑗 are the liquid viscosity, 𝑙 𝑖𝑗 pore throat length, and pore throat radius, respectively. 6 Gas transport: The transport of gas through pore throats is complicated: because the partial pressure of the solvent vapor is comparable to that of the purge gas, a multicomponent transport model is required. Here, we adopt the dusty gas model (DGM), which have been widely used for multicomponent mass transfer in nanopores. The transport of an ideal mixture of d gas species through a pore with a radius of 𝑟𝑖𝑗 can be described using30-31 ∑𝑑𝑛=1 𝑛≠𝑚 𝑥𝑛 𝐽𝑚 −𝑥𝑚 𝐽𝑛 +𝜌 𝜌𝑡 𝐷𝑚𝑛 𝐽𝑚 𝑡 𝐷𝑚𝑘 1 𝑑𝑝𝑚 = −𝑝 𝑑𝑧 −𝜇 𝑥𝑚 𝑘𝑝 𝑑𝑝 𝑚𝑖𝑥 𝐷𝑚𝑘 𝑑𝑧 (𝑚 = 1,2, … 𝑑) (7) where 𝐽𝑚 , 𝜌𝑚 , and 𝑥𝑚 are the molar flux, averaged molar density, and molar ratio of specie m, respectively. 𝑝 is the pressure of the 𝑑-component gas mixture and 𝑝𝑚 is the partial pressure of a specie m. 𝜌𝑡 is the total molar density, 𝑘𝑝 is the effective permeability, 𝜇𝑚𝑖𝑥 is the mixture viscosity, 𝐷𝑚𝑛 is the mutual diffusivity between species m and n. 𝐷𝑚𝑘 = 2𝑟𝑖𝑗 3 8𝑅𝑔 𝑇 √ 𝜋𝑀 is the Knudsen diffusivity of species m (𝑀𝑚 is the molar mass, 𝑚 Rg is the ideal gas constant and 𝑇 is the temperature). When applied to a binary mixture of purge gas (species 1) and solvent vapor (species 2) that behaves ideally, the molar flux of each gas species can be reorganized and simplified as 𝐽1 = 𝜌1 𝑢1 = − (( 𝐽2 = 𝜌2 𝑢2 = − ((𝑥 (𝑥1 𝐷2𝑘 𝐷1𝑘 +𝐷12 𝐷1𝑘 ) ⁄𝑅𝑔 𝑇 𝑥2 𝐷1𝑘 +𝑥1 𝐷2𝑘 +𝐷12 𝑥2 𝐷2𝑘 𝐷1𝑘 /𝑅𝑔 𝑇 2 𝐷1𝑘 +𝑥1 𝐷2𝑘 +𝐷12 𝑘𝑝 𝜌2 +𝜇 𝑚𝑖𝑥 ) 𝑑𝑝1 𝑑𝑧 𝑘𝑝 𝜌1 +𝜇 +( 𝑚𝑖𝑥 ) 𝑑𝑝1 𝑑𝑧 + (𝑥 𝑥1 𝐷2𝑘 𝐷1𝑘 /𝑅𝑔 𝑇 2 𝐷1𝑘 +𝑥1 𝐷2𝑘 +𝐷12 (𝑥2 𝐷2𝑘 𝐷1𝑘 +𝐷12 𝐷2𝑘 ) ⁄𝑅𝑔 𝑇 𝑥2 𝐷1𝑘 +𝑥1 𝐷2𝑘 +𝐷12 The overall gas volume flux 𝑄𝑖𝑗𝑔 through pore throat ij can be wrote as 𝑔 𝑄𝑖𝑗 = (𝐷2𝑘,𝑖𝑗 𝐷1𝑘,𝑖𝑗 +𝐷12,𝑖𝑗 𝐷2𝑘,𝑖𝑗 )⁄ 𝑅𝑔 𝑇 where 𝐾𝑖𝑗𝐴 = ( 𝑥2,𝑖𝑗 𝐷1𝑘,𝑖𝑗 +𝑥1,𝑖𝑗 𝐷2𝑘,𝑖𝑗 +𝐷12,𝑖𝑗 2 𝜋𝑟𝑖𝑗 (𝐽1,𝑖𝑗 +𝐽2,𝑖𝑗 ) 𝑔 𝑘𝑝 𝜌𝑖𝑗 +𝜇 𝑚𝑖𝑥 𝑔 ) 𝜌𝑖𝑗 2 𝜋𝑟𝑖𝑗 𝑔 𝜌𝑖𝑗 𝑙𝑖𝑗 𝑘𝑝 𝜌2 +𝜇 𝑚𝑖𝑥 ) 𝑘𝑝 𝜌1 +𝜇 𝑑𝑝2 𝑑𝑧 𝑚𝑖𝑥 ) 𝑑𝑝2 𝑑𝑧 ) ) (8b) 𝑔2 𝑔2 𝑔 𝑔 = − (𝐾𝑖𝑗𝐴 (𝑃𝑖 − 𝑃𝑗 ) + 𝐾𝑖𝑗𝐵 (𝑃𝑖 − 𝑃𝑗 )) and 𝐾𝑖𝑗𝐵 = (𝑥 𝐷12,𝑖𝑗(𝐷2𝑘,𝑖𝑗 −𝐷1𝑘,𝑖𝑗 ) ⁄𝑅𝑔 𝑇 2,𝑖𝑗 𝐷1𝑘,𝑖𝑗 +𝑥1,𝑖𝑗 𝐷2𝑘,𝑖𝑗 +𝐷12,𝑖𝑗 ) 2 𝜋𝑟𝑖𝑗 𝑔 𝜌𝑖𝑗 𝑙𝑖𝑗 (8a) (9) are transport coefficients. 𝑃𝑖𝑔2 is the vapor’s partial pressure in pore body i, 𝜌𝑖𝑗𝑔 = 0.5(𝜌𝑖𝑔 + 𝜌𝑗𝑔 ) is the average molar density in pores i and j. Substituting Eqs. (6) and (9) into Eq. (5), we obtain 𝑔2 𝑔2 𝑔2 𝑔 𝑔 𝑔 𝑖 ∑𝑁 (𝐾𝑖𝑗𝑙 (𝑃𝑖𝑙 − 𝑃𝑗𝑙 ) + 𝐾𝑖𝑗 (𝑃𝑖 − 𝑃𝑗 ) + 𝐾𝑖𝑗 (𝑃𝑖 − 𝑃𝑗 )) = 0 𝑗=1 Eq. (10) can be reformulated using the average pressure of the pore body 𝑃̃𝑖 = (1 − 𝑆𝑖𝑙 )𝑃𝑖𝑔 + 𝑆𝑖𝑙 𝑃𝑖𝑙 as 𝑁𝑖 𝑖 ∑𝑁 (𝐾𝑖𝑐 𝑃𝑖𝑐 − 𝐾𝑗𝑐 𝑃𝑗𝑐 + 𝐾𝑖𝑗𝑔2 (𝑃𝑖𝑔2 − 𝑃𝑗𝑔2 )) (𝐾𝑖𝑗𝑙 + 𝐾𝑖𝑗𝑔 )(𝑃̃𝑖 − 𝑃̃𝑗 ) = − ∑𝑗=1 𝑗=1 where 𝐾𝑖𝑐 = 𝐾𝑖𝑗𝑔 𝑆𝑖𝑙 − 𝐾𝑖𝑗𝑙 (1 − 𝑆𝑖𝑙 ) and 𝐾𝑗𝑐 = 𝐾𝑖𝑗𝑔 𝑆𝑗𝑙 − 𝐾𝑖𝑗𝑙 (1 − 𝑆𝑗𝑙 ). 7 (10) (11) Using the kinetic theory of gas transport inside nanopores and considering the first order slippage effect, the effective permeability for gas through a pore throat ij is given by32 𝑘𝑝 = 2 𝑟𝑖𝑗 8 2 (1 + 4Kn(𝜎 − 1)) 𝑣 (12) where 𝜎𝑣 is the tangential momentum accommodation coefficient. We assume 𝜎𝑣 = 1, which corresponds to a rough pore surface that reflects all molecules diffusively. Kn = λ/2𝑟𝑖𝑗 is the Knudsen number of the gas mixture (λ is the mean free path for the gas mixture, see below), which defines the ratio of the Knudsen diffusion to self-diffusion in a bulk gas. Properties of gas mixture: To complete the above gas transport model, the mutual diffusivity, mixture viscosity, and the mean free path of the gas mixture are needed. The mutual diffusivity in a pore throat ij is approximated using the empirical formula30, 33 𝑥1 𝐷12 = (𝐷𝑥1→1 ) (𝐷𝑥1→0 ) (1−𝑥1 ) (13) where 𝑥1 = 0.5(𝑥1,𝑖 + 𝑥1,𝑗 ) is the molar fraction of species 1 in a pore throat ij, which is taken as the average of the molar faction of species 1 in pore body i (𝑥1,𝑖 ) and pore body j (𝑥1,𝑗 ). The bracketed terms are the infinite dilution values for the Maxwell-Stefan diffusivity at either end of the composition range and are 2 √𝑚𝑖𝑘𝐵 𝑇 ⁄𝜋 given by 𝐷𝑥𝑖→1 = 3 𝑀𝑖𝜎𝑎𝑖 𝜌𝑖 (𝑖 = 1,2), where 𝑀𝑖 and 𝜎𝑎𝑖 are the molar mass and effective collision cross- section area of specie i, respectively. We approximate the viscosity of the binary gas mixture using an empirical correlation:34 𝜇𝑚𝑖𝑥 = 𝜇1 𝑥 1.385𝜇1 1+ 2 𝑥1 𝐷12 𝜌1 + 𝜇2 𝑥 1.385𝜇2 1+ 1 where 𝜇𝑖 is the viscosity for a pure species i. The mean free path of a binary gas mixture is given by 𝜆= 𝑥 𝑥1 𝑚1+𝑥2 𝑚2 𝑥 (14) 𝑥2 𝐷12 𝜌2 √2𝜎𝑎11 𝜎𝑎22 (𝑀1 𝜌1 +𝑀2 𝜌2 ) (15) where 𝑚𝑖 (𝑖 = 1,2) are molecular mass for the species i. B.3 Liquid vaporization model The evaporation of liquid solvents at the corner of pore bodies and the vapor’s subsequent transport to the pore center are essential steps of the vaporization process induced by purge gas. Since we assume that drying proceeds under the isothermal condition, the vaporization rate is limited by the vapor transport from the liquid surface toward the pore’s center rather than the kinetics of liquid evaporation. Therefore, the 8 vapor density on the liquid surface is equal to the liquid’s saturation vapor density 𝜌𝑠𝑎𝑡 and the vaporization from the liquid surface inside a pore body i is given as: 𝑒𝑣𝑝 𝑚̇𝑖 𝑔2 = 𝛽𝐴𝑠 (𝜌𝑠𝑎𝑡 − 𝜌𝑖 ) (16) where 𝛽 is the mass transfer coefficient. 𝐴𝑠 is the interfacial area between the liquid solvent and the gas phase inside the pore body and it is given by:24 𝐴𝑠 = 2 144 2 𝑟 (1 − 𝑆𝑖𝑙 )3 𝜋 𝑖 { 12𝜋𝑟𝑖 𝛾 8𝜋𝛾2 𝑃𝑐 − , 𝑃𝑐2 , 𝑆𝑖𝑙 ≥ 0.476 (17) 𝑆𝑖𝑙 < 0.476 The mass transfer of vapor from the surface of liquids at pore corners to the pore interior is characterized using the Sherwood number 𝑆ℎ = 2𝛽𝑟𝑖 /𝐷𝑠, which depends on the pore’s shape. Here 𝑆ℎ is taken as 2.98, which corresponds to square-shaped pores.35 With the above vaporization model, the mass balance of the solvent vapor in a pore body i is given by 𝑔 𝑁 𝑔2 𝑔2 𝑔2 𝑖 𝑉𝑖 Δ𝜌𝑖 /∆𝑡 = ∑𝑗=1 𝜋𝑟𝑖𝑗2 𝐽𝑖𝑗 + 𝐴𝑠 𝛽(𝜌𝑠𝑎𝑡 − 𝜌𝑖 ) (18a) 𝑔2 𝑔2 𝑔2 𝑔 𝑔 𝑔2 𝑔 𝐽𝑖𝑗 = − (𝐶𝑖𝑗 (𝑃𝑖 − 𝑃𝑗 ) + 𝐶𝑖𝑗 (𝜌𝑖 − 𝜌𝑗 )) (18b) where 𝑉𝑖 𝑔 = 𝑉𝑖 (1 − 𝑆𝑖𝑙 ) is the volume of the gas phase inside a pore body i and ∆𝑡 is the time step (see B.5 for details), 𝐶𝑖𝑗𝑔 = 𝑥 coefficients. 𝑥2,𝑖𝑗𝐷1𝑘,𝑖𝑗 𝐷2𝑘,𝑖𝑗 /𝑅𝑔 𝑇 1,𝑖𝑗 𝐷2𝑘,𝑖𝑗 +𝑥2,𝑖𝑗 𝐷1𝑘,𝑖𝑗 +𝐷12,𝑖𝑗 𝑘𝑝 +𝜇 𝑚𝑖𝑥 𝑔2 and 𝐶𝑖𝑗𝑔2 = 𝑥 𝜌𝑖𝑗 𝐷2𝑘,𝑖𝑗 𝐷12,𝑖𝑗 1,𝑖𝑗 𝐷2𝑘,𝑖𝑗 +𝑥2,𝑖𝑗 𝐷1𝑘,𝑖𝑗 +𝐷12,𝑖𝑗 are transport B.4 Initial and boundary conditions Initially, the pore network is fully saturated with the liquid solvents: 𝑆𝑖𝑙 (𝑡 = 0) = 1.0 (19) During drying, the pressure at the pore network’s inlet (y=0 in Fig. 2) is fixed to that of the upstream purge gas reservoir 𝑃𝑢 and no solvent vapor is transported into the pore from the upstream reservoir: 𝑦 = 0: 𝑃𝑔 |𝑦=0 = 𝑃𝑢 (20a) 𝑦 = 0− : 𝐽𝑔2 = 0 (20b) 𝑦 = 𝐿: 𝑃𝑔 |𝑦=𝐿 = 𝑃𝑑 (20c) The pressure at the pore network’s outlet (𝑦 = 𝐿 in Fig. 2) is fixed to that at the downstream reservoir (𝑃𝑑 ): 𝑦 = 𝐿: 𝜕𝜌 𝑔2 𝜕𝑦 =0 (20d) 9 B.5 Rules 1. Threshold pressure for pore throats. Because the contact angle of solvent liquids on the solid surfaces is zero, the threshold pressure for pore throats filled with liquids is 2𝛾 𝑃𝑡ℎ,𝑖𝑗 = 𝑟 (21) 𝑖𝑗 A pore throat filled with liquid can be invaded by the gas phase if the capillary pressure in the neighboring pore body exceeds the threshold pressure of this pore throat: 24 𝑃𝑐 (𝑆𝑖𝑙 ) > 𝑃𝑡ℎ,𝑖𝑗 or 𝑃𝑐 (𝑆𝑗𝑙 ) > 𝑃𝑡ℎ,𝑖𝑗 (22) Because the liquid saturation 𝑆𝑖𝑙 in a pore body can be reduced by vaporization of the liquid solvents in it, Eq. (22) implies that evaporation can indirectly affect the invasion of a pore throat by liquid solvents or gas. 2. Re-imbibition of pore throats. If a pore body i is refilled with liquid solvent (𝑆𝑖𝑙 = 1.0), an empty throat (i.e., a throat free of liquid) ij can be re-imbibed with the liquid under the condition29 𝑔 𝑃𝑖 > 𝑃𝑗 𝑔 (23) where 𝑃𝑖𝑔 is the gas pressure of the target pore body and 𝑃𝑗𝑔 is the gas pressure in the pore body connected to the target pore body by the throat ij. 3. Selection of time step. The time step for each iteration in the pore network simulation is determined based on the local time steps required for draining a pore body to the critical state, i.e., when the local capillary pressure satisfies 𝑃𝑐 (𝑆𝑖𝑙 ) = min{𝑃𝑡ℎ,𝑖𝑗 }. Since the imbibition process is allowed to occur locally in pore bodies, the local time step can be determined by the time required to refill the pore body 24: 𝑉𝑖 𝑁 𝑙 ∑ 𝑖 𝑄𝑖𝑗 𝑗=1 ∆𝑡𝑖 = { 𝑉𝑖 𝑁 𝑙 ∑ 𝑖 𝑄𝑖𝑗 𝑗=1 𝑁𝑖 𝑙 (𝑆𝑖𝑙 − 𝑆𝑖,𝑐𝑡 ) local drainage, ∑𝑗=1 𝑄𝑖𝑗 <0 (1 − 𝑆𝑖𝑙 ) 𝑁𝑖 𝑙 local imbibition, ∑𝑗=1 𝑄𝑖𝑗 >0 (24) where 𝑆𝑖,𝑐𝑡 is the saturation corresponding to the critical state. Using the 𝑃𝑐 − 𝑆𝑖𝑙 relation by Eq.(4), the critical saturation can be given as 1 𝑆𝑖,𝑐𝑡 = − 6.83 ln(1 − 𝑟 2𝛾 ) 𝑖 min{∆𝑃𝑢𝑑 ,min{𝑃𝑡ℎ,𝑖𝑗 }} (25) where ∆𝑃𝑢𝑑 = 𝑃𝑢 − 𝑃𝑑 is the global pressure difference between the upstream and downstream reservoir. min{𝑃𝑡ℎ,𝑖𝑗 } is the minimum threshold pressure among the liquid throats connected to the target pore body. Then the global time step is selected as the minimum of all local time steps, i.e., ∆𝑡 = min{∆𝑡𝑖 }. 10 C. Computational algorithms Starting from the initial conditions given by Eq. (19), Eqs. (4) and (11) are combined with the boundary conditions given by Eq. (20) to solve the pressure field in the pore network using the direct sparse-matrix solver (DSS).36 Using this pressure field, the solvent vapor density for the next time step is computed by solving Eqs. (17) and (18) subject to the boundary conditions in Eq. (20) using an implicit scheme. These calculations are repeated till the entire pore network is free of liquid solvents. These and other details of the computational algorithms for solving the pore network model are summarized in Fig. 3. Figure 3. The computational procedure and algorithm for solving the pore network model in which drainage and evaporation can occur concurrently. III. Results and Discussion We consider the drying of a porous coal cake initially saturated with liquid solvent pentane by the flow of a purge gas N2 through the cake. Because coal particles’ diameter is often less than a few micrometers, the size of the pores in the coal cake can be a few hundred nanometers. Here, the coal cake is represented using a pore network consisting of 80 × 80 pore bodies, with each pore body connected to its 4 nearest pore bodies as shown in Fig. 2a. The size of pore bodies follows the truncated log-normal distribution with 𝜎𝑏 = 0.22, 𝑟𝑏,𝑚 = 1.2 μm, 𝑟𝑏,𝑚𝑖𝑛 = 0.8 μm and 𝑟𝑏,𝑚𝑎𝑥 = 1.4 μm . The size of pore throats follows the same distribution with 𝜎𝑡 = 0.18 , 𝑟𝑡,𝑚 = 300 nm, 𝑟𝑡,𝑚𝑖𝑛 = 200 nm and 𝑟𝑡,𝑚𝑎𝑥 = 400 nm . The size distributions of 11 pore bodies and throats are shown in Fig. 2c and 2d. The pore bodies and throats generated above are randomly packed into the pore network. The system temperature is 300 K. The material properties of the purge gas (N2) and solvent (pentane) at this temperature are summarized in Table 1. Note that the effective cross-section areas of the N2 and pentane molecules (a,1 and a,2) are calculated based on the results in Ref. 37. Table 1. Properties of the purge gas N2 (species 1) and solvent (species 2). m1 = 4.65×10-26 kg a,1 = 0.43 nm2 m2 =12.0×10-26 kg a,2 = 2.09 nm2 sat,2 = 28.1 mol/m3 1 = 17.8 Pa∙s (gas state) Psat,2 = 0.7 bar 2 = 70.0 Pa∙s (gas state) 2= 626.0 kg/m (liquid state) 3 1 = 224.0 Pa∙s (liquid state)  = 15 mN/m The pressure in the reservoir downstream the pore network, Pd, is 1.0 bar. For selection of the upstream gas pressure (Pu) , we identify a minimal pressure difference as 𝑚𝑖𝑛 ∆𝑃𝑢𝑑 =𝑟 2𝛾 𝑡,𝑚𝑖𝑛 (28) When 𝑃𝑢 − 𝑃𝑑 = ∆𝑃𝑚𝑖𝑛 𝑢𝑑 , not all all pores inside the pore network can be emptied by drainage because the pressure drop along individual throats is smaller than the global pressure difference across the entire pore 𝑚𝑖𝑛 network. Nevertheless, ∆𝑃𝑢𝑑 provides an indication of the ability of the applied pressure difference to drain 𝑚𝑖𝑛 liquids from the pore network. For the pore network studied here, ∆𝑃𝑢𝑑 =1.5 bar. In our simulations, we select Δ𝑃𝑢𝑑 = 𝑃𝑢 − 𝑃𝑑 = 1.4, 1.5, and 2.0 bar to study how the applied pressure difference affects the drying behavior. These values are within the range used in the dewatering-by-displacement technology.22 A. Macroscopic drying behaviors To quantify how a porous medium initially saturated with liquid solvents is dried, we define a degree of drying 𝛿 using 𝛿(𝑡) = 1 − 𝑆̅(𝑡), where 𝑆̅ is the ratio of the liquid solvent mass at a time 𝑡 and the initial solvent mass in the porous medium. In purge gas-assisted drying of a porous medium, solvent can be removed as liquid through drainage (termed drainage effect) and by the transport of vaporized solvents out ̇ of the porous medium (termed evaporation effect). We thus decompose the net solvent removal rate 𝐹𝑛𝑒𝑡 ̇ = 𝐹𝑑𝑟𝑛 ̇ + 𝐹𝑒𝑣𝑝 ̇ , where 𝐹𝑑𝑟𝑛 ̇ ̇ are solvent fluxes at the porous medium’ from the porous medium as 𝐹𝑛𝑒𝑡 and 𝐹𝑒𝑣𝑝 outlet due to the drainage and evaporation effects, respectively. To quantify the relative importance of the drainage and evaporation effects, we further define an evaporation-to-drainage ratio as 𝛼 = 𝐹̇ 𝑒𝑣𝑝 /𝐹̇ 𝑑𝑟𝑛 . 12 Figure 4. The overall drying behavior at different applied pressure Δ𝑃𝑢𝑑 . (a-b) The evolution of the ratio of removal rates of solvent as vapor and liquid (a) and the total solvent removal rate (b) as a function of the degree of drying . (c) The evolution of the average liquid saturation in the porous medium. The dashed line corresponds to the simulation in which evaporation is disabled. Figure 4a shows the evolution of 𝛼 during drying under three pressure differences Δ𝑃𝑢𝑑 . As drying proceeds, 𝛼 generally increases, i.e., the evaporation effect becomes more and more important. Base on the relative importance of the drainage and evaporation effects, the drying process under all Δ𝑃𝑢𝑑 can be divided into three stages: stage I with 𝛼 < 0.3, during which the solvent removal is dominated by the drainage effect; stage II with 0.3 < 𝛼 < 3.0 , during which the drainage and evaporation effects are comparable; stage III with 𝛼 mostly larger than 3, during which the drying is largely dominated by the evaporation effects. For all Δ𝑃𝑢𝑑 , the degree of drying spanned by each stage is similar, i.e., 𝛿 ≲ 0.65 in stage I, 0.65 ≲ 𝛿 ≲ 0.82 in stage II, and 𝛿 ≳ 0.82 in stage III. During stage III, 𝛼 fluctuates notably and can reach ~0.3-1 (i.e., drainage notably contributes to solvent removal), and this is especially pronounced for the smallest Δ𝑃𝑢𝑑 studied (1.4 bar). The contribution of drainage to the solvent removal during each stage is summarized in Table 2. We observe that, in the first two stages, the contribution of drainage to solvent removal is similar for all Δ𝑃𝑢𝑑 . However, in stage III, drainage contributes more significantly at lower Δ𝑃𝑢𝑑 , e.g., the contribution of drainage to solvent removal increases from 6.2% at Δ𝑃𝑢𝑑 = 2.0 bar to 21.4% at Δ𝑃𝑢𝑑 = 1.4 bar. Nevertheless, because most of the solvents are removed during stage I, the overall contribution of drainage to solvent removal does not differ greatly for the Δ𝑃𝑢𝑑 studied: as Δ𝑃𝑢𝑑 decreases from 2.0 to 1.4 bar, the contribution of drainage increases only by 1.5%. Nevertheless, because solvent removal by drainage avoids the energy cost due to the heat of vaporization needed in evaporation-induced solvent removal, even a small increase of drainage’s contribution is desirable, especially in large-scale applications. ̇ during drying. Overall, 𝐹𝑛𝑒𝑡 ̇ is higher Figure 4b shows the evolution of the net solvent removal rate 𝐹𝑛𝑒𝑡 ̇ curves at different Δ𝑃𝑢𝑑 exhibit similarities. First, 𝐹𝑛𝑒𝑡 ̇ fluctuates notably except at for larger Δ𝑃𝑢𝑑 . The 𝐹𝑛𝑒𝑡 ̇ is similar: in stage I, 𝐹𝑛𝑒𝑡 ̇ increases first and then the beginning of stage I. Second, the time evolution of 𝐹𝑛𝑒𝑡 ̇ is ̇ continues to decrease in stage II but plateaus in stage III. Significant fluctuation of 𝐹𝑛𝑒𝑡 decreases. 𝐹𝑛𝑒𝑡 observed in the stage III, especially at Δ𝑃𝑢𝑑 = 1.4 bar. 13 Table 2. The contribution of drainage to the solvent removal and the duration of each drying stage. Δ𝑃𝑢𝑑 = 1.4 𝑏𝑎𝑟 Δ𝑃𝑢𝑑 = 1.5 𝑏𝑎𝑟 Δ𝑃𝑢𝑑 = 2.0 𝑏𝑎𝑟 stage II 57.8% 53.1% 61.0% stage III 21.4% 8.9% 6.2% total 75.9% 73.3% 74.4% stage I 46.82 39.58 27.77 stage II 53.91 47.77 32.22 stage III 139.30 110.50 68.46 total 240.03 197.85 128.45 stage I contribution of drainage to the solvent removal duration (ms) 96.1% 96.6% 96.9% Figure 4c shows the time evolution of the averaged solvent liquid saturation in the porous medium. The liquid saturation decreases sharply during stage I, but the decrease of saturation slows down in stage II and even more so in stage III. These observations are consistent with the fact that, as drying proceeds, drying becomes controlled more by the evaporation effects and the solvent removal rate decreases (see Fig. 4a and 4b). Note that the evaporation effects are essential for the complete removal of solvents. To appreciate this, we perform simulations in which the vaporization of liquid solvents is disabled and only drainage is allowed. Comparison of the result of this simulation (the dashed line in Fig. 4c) with the above result shows that, as the average liquid saturation reduces to ~0.3, the solvent removal starts to be affected by the evaporation effect. In absence of the evaporation effects, a significiant fraction of the liquid solvents remains trapped in the porous medium. Figure 4c shows that the evaporation effects-dominated stage III lasts longer compared to the other two stages, i.e., the throughput of drying is limited by the removal of solvents by evaporation effects. This is seen more clearly in Table 2, in which the time corresponding to each drying stage is listed. The time needed for the stage III is especially long under the applied pressure difference of 1.4 bar, which is partly caused by the smaller driving force for gas transport through the pore network. B. Microscopic drying processes The existence of three drying stages, the contributions and interplay of drainage and evaporation effects during drying, and finally the evolution of the drying rate revealed in Fig. 4 can be understood by examining the microscopic processes in the pore network during drying. Figure 5 shows the snapshots of the liquid saturation distribution in the pore network at four representative degrees of drying (𝛿 ) when the applied pressure difference is ∆𝑃𝑢𝑑 = 1.5 bar. Similar results are observed for other ∆𝑃𝑢𝑑 and are not shown here. 14 Once the pressure difference is imposed across the pore network, the purge gas starts to invade into the pore network. The purge gas travels preferentially through pathways with wider pore throats and displaces the liquids in a piecemeal manner. This is evident in Fig. 5a ( 𝛿 =0.16), where finger-shaped gas paths are observed. The formation of these paths, often referred to as viscous fingering, 38-39 is generally considered as the onset and evolution of instabilities when a more viscous fluid (here, liquid solvent) is displaced by a less viscous fluid (here, purge gas). As the viscous fingers grow, the network formed by liquid-saturated pore bodies is fragmented into liquid clusters (hereafter, a liquid cluster is defined as a collection of pore bodies fully saturated by liquids and connected continuously by liquid-filled throats). In particular, the initial main liquid cluster spanning the entire pore network is fragmented into small liquid clusters, which are broken into even smaller clusters later. Many small liquid clusters, along with some large liquid clusters, appear inside the system (see Fig. 5b, where 𝛿 =0.5). During the early part of this process, as viscous fingers move toward but have not yet reached the pore network’s outlet, the resistance for drainage decreases and the drainage rate increases, which is consistent with the initial increase of the drying rate (see Fig. 4b, 𝛿 ≲ 0.2 − 0.3). After the purge gas breaks through the pore network, drainage from the pore network’s outlet becomes more limited. Therefore, the drying rate decreases as shown in Fig. 4b and the evaporation effects begin to contribute to solvent removal. However, until drainage pathways are diminished by the fingering of gas pathways, the removal of solvents is dominated by drainage and drying is in the stage I identified in Fig. 4a. Figure 5. The evolution of liquid saturation in the pore network when the applied pressure difference Δ𝑃𝑢𝑑 is 1.5 bar. As drying proceeds, the liquid saturation throughout the pore network reduces, and this is evident in the snapshot in Fig. 5c (𝛿 =0.69). Consequently, the drainage of liquid solvents through the network is reduced while more pathways for vaporized solvents to be “flushed” out of the network emerge. The solvent removal by the evaporation effects thus increases. The contributions of drainage and evaporation effects eventually becomes comparable and drying enters stage II. Because drainage is much more effective in removing solvent (note that the density of liquid solvents is >100 times larger than the density of vaporized solvents), the net drying rate continues to decrease in stage II. 15 In stage II of the drying process, a drying front, behind which few liquid clusters exist, appears near the pore network’s inlet (see Fig. 5c) and moves downstream as drying proceeds. However, it important to note that solvents are removed from individual pores both near and far ahead of this drying front. As drying continues, more and more pores near the network’s outlet become only partially occupied by liquid solvents. Therefore, solvent removal by drainage (evaporation effects) diminishes (becomes significant), and drying enters stage III as shown in Fig. 4a. In stage III, the drying rate is quite stable for Δ𝑃𝑢𝑑 = 1.5 and 2.0 bar (see Fig. 4b). Even for Δ𝑃𝑢𝑑 = 1.4 bar, the drying rate largely fluctuates around a constant value. These trends can be understood as follows. At this stage, small liquid clusters are mostly trapped in pores with narrow throats. The pathway for gas (purge gas + solvent vapor) transport evolves only slowly and the gas flow rate is relatively stable. Because the gas phase at the pore network’s exit is generally saturated with the solvent vapor (except toward the very end of drying), the net flux of vaporized solvents out of the porous medium is rather stable. Because drying is dominated by the vaporized solvents, the total drying rate is relatively stable. C. The interplay between drainage and evaporation effects The above analysis highlights the role of drainage and evaporation effects in purge gas-assisted drying of porous media and helps understand the main features of the macroscopic drying behavior. However, some features of the drying data shown in Fig. 4 still cannot be readily understood. For example, a notable drainage flux can still appear after the evaporation effects start to dominate drying (see Fig. 4a and Table 2), which is also manifested as the spikes in the net drying rate in stage III (see Fig. 4b). These features essentially arise from the coupling between the drainage and evaporation effects. To understand the coupling between drainage and evaporation effects and its impact on the macroscopic drying behavior, it is instructive to examine the evolution of liquid clusters in the pore network, which are affected oppositely by the two processes: in drainage, purge gas breaks through liquid clusters to fragment them into small clusters; the evaporation in pore bodies tends to eliminate small liquid clusters. Figure 6a shows the evolution of the number of liquid clusters 𝑁𝑐 in the pore network. In stage I, evaporation effects are minor and 𝑁𝑐 increases rapidly as viscous fingers break through increasingly smaller clusters. In stage II, fragmentation of liquid clusters by purge gas “fingers” becomes weaker while the elimination of liquid clusters by evaporation effects becomes significant. The competition of these processes leads to a plateau of 𝑁𝑐. In stage III, where evaporation effects dominate, 𝑁𝑐 decreases. 16 Figure 6. (a) The evolution of the number of liquid clusters in the pore network during drying and drainage-only operations. (b-c) The distribution of the liquid saturation in a pore network at a degree of drying of 𝛿 = 0.7 as obtained from drying simulations (b) and drainage-only simulations (c). The applied pressure difference ∆𝑃𝑢𝑑 is 1.4 bar. To more directly appreciate the coupling between drainage and evaporation effects, we performed a new simulation in which the evaporation effects are turned off so that only drainage exists. The setup of the pore network is identical to the drying simulation reported above and the pressure difference Δ𝑃𝑢𝑑 is set to 1.4 bar. Figure 6a shows that, in absence of evaporation effects, the evolution of 𝑁𝑐 is similar in the drainage and drying simulations till the degree of drying 𝛿 reaches 0.6, which is expected because evaporation effects are weak in the stage I drying. However, in absence of evaporation effects, 𝑁𝑐 increases with 𝛿 till it reaches ~0.7. After that, solvents are permanently trapped in small liquid clusters (see Fig. 6b) and 𝑁𝑐 remains constant. Importantly, in this case, numerous liquid clusters with small size are scattered throughout the pore network (see Fig. 6b), in sharp contrast with the fact that, at the same 𝛿 , the liquid clusters in the simulation with both drainage and evaporation effects have much larger size and are distributed closer to the pore network’s outlet (see Fig. 6c). This difference can be understood as follows. Once evaporation enables the emptying of a pore throat connected to a liquid cluster, drainage can displace the liquid in the cluster to neighboring pore bodies. This effectively drives the liquid cluster downstream and even allows it to merge with other liquid clusters and form the large liquid clusters seen in Fig. 6c. The coupling between drainage and evaporation effects revealed through the comparison between Fig. 6b and 6c affects the macroscopic drying behavior. For example, the liquid clusters near the pore network’s outlet can be removed through drainage, which helps explain why the solvent flux at the porous medium’s outlet frequently includes a significant faction of liquid even in stage III of drying (see Fig. 4a, especially when Δ𝑃𝑢𝑑 is 1.4 bar). The extent of the coupling between drainage and evaporation effects during a drying process is affected by both the operating conditions (e.g., the applied pressure difference) and the stage of drying. To gauge the extent of this coupling, we note that, the coupling between the drainage and evaporation effects can be 17 considered as three steps looped together. First, evaporation triggers the emptying of pore throats (see Eq. (22)) and prompts liquid to drain from pores. Second, drainage generates new gas transport pathways and affects the liquid saturation distribution in the pore network. Third, the altered gas transport pathways and liquid saturation distribution in turn affect the evaporation of liquid solvents in the partially saturated pores. These looped steps are most easily delineated for a liquid cluster surrounded by very narrow pore throats (see Appendix) but they are applicable to all liquid clusters in a pore network. Because the coupling loop between drainage and evaporation effects is triggered when evaporation induces the emptying of pore throats, we can infer the extent of such coupling by counting the number of its triggering events (i.e., evaporation-induced emptying of pore throat) 𝑁𝑡𝑟𝑖𝑔 . A larger 𝑁𝑡𝑟𝑖𝑔 corresponds to a more extensive coupling between the drainage and evaporation effects during drying. Figure 7. The statistics of the triggering events (evaporation-induced emptying of pore throats) in the drainageevaporation coupling loops during drying. (a-c) The histogram of the number of triggering events as a function of the degree of drying 𝛿 under different applied pressure differences. Figure 7a, 7b, and 7c show the histogram of 𝑁trig under an applied pressure difference of ∆𝑃𝑢𝑑 = 1.4, 1.5, and 2.0 bar, respectively. There are few triggering events for drainage-evaporation coupling in stage I, but these events start to appear in stage II, when the drainage and evaporation effects become comparable. In stage III, when the evaporation effect dominates, many triggering events are detected. The number of triggering events is largest at ∆𝑃𝑢𝑑 = 1.4 bar, i.e., the coupling between the drainage and evaporation effects is most extensive at this applied pressure difference. As discussed in section III.B, liquid clusters can be removed in two processes in presence of the coupling effects: (1) liquid clusters can migrate downstream, merge with other liquid clusters, and be drained out of the pore network at its outlet; (2) liquid clusters can be directly removed from the porous media via evaporation effects. Because the coupling of drainage and evaporation effects is more extensive for ∆𝑃𝑢𝑑 = 1.4 bar, the first process is strongest at this applied pressure difference. This explains why, during stage III of drying, the removal of liquid solvents by drainage from the pore network’s outlet is more significant for ∆𝑃𝑢𝑑 = 1.4 bar than for ∆𝑃𝑢𝑑 = 1.5 and 2.0 bar (see Fig. 4a and Table 2). 18 IV. Conclusions In this work, we develop a new pore network model for purge gas-assisted drying of nanoporous media. The model considers multiphase and multiphysics processes such as liquid drainage, evaporation, and transport of gas mixtures through nanopores. Solutions of the model indicate that the drying process can be divided into three stages. In stage I, drying is dominated by drainage, in which gas flow displaces the liquid solvent in a fingering pattern. The net drying flux increases before the purge gas breaks through the porous medium and decreases after that. During this process the number of liquid clusters in the porous medium increases. In stage II, the evaporation effects become comparable to the drainage effect. While drainage tends to increase the number of liquid clusters in the porous medium, evaporation tends to reduce the number of liquid clusters. Hence, the number of liquid clusters in the porous medium plateaus. As drying proceeds to stage III, the liquid removal is largely dominated by the evaporation effects, although notable drainage can also occur at the porous medium’s outlet, especially when the pressure driving the purge gas is low. A key feature of the present drying process is the coupling between the drainage and evaporation effects. Evaporation can trigger three-step coupling loops that begin with the emptying of pore throats, continue with the drainage of liquids and fragmentation of liquid clusters, and continue with evaporation in the newly created partially saturated pore bodies. Because of the interplay between drainage and evaporation effects in these coupling loops, studying purge gas-assisted drying by dividing it into a drainage period and an evaporation period and studying them separately is generally inadequate (cf. the vastly different liquid saturation distributions in Fig. 6b and 6c). The coupling between drainage and evaporation effects depends on operating conditions, e.g., it is extensive in a porous medium when the applied pressure difference is comparable to the minimum threshold pressure for gas to invade the narrowest pore throats. The extent of coupling greatly affects the evolution of the liquid clusters in the porous medium and consequently the drying behavior, e.g., as the coupling becomes more extensive, liquid removal through drainage in stage III of the drying process increases, which helps increases the energy efficiency of drying. The insight on the coupling between the drainage and evaporation effects in purge gas-assisted drying helps guide future application of this method. Appendix In the main text, the coupling of drainage and evaporation effects is conceptualized as three looped steps occurring in a space initially occupied by a liquid cluster: evaporation-induced pore throat emptying, drainage through the emptied throats and alteration of liquid saturation distribution, and alteration of the evaporation in that space. Here we perform simulations to illustrate these steps in a pore network where a 19 well-defined liquid cluster exists. The statistics of the pore bodies and throats of this pore network are the same as those in the main text except that the radius of the pore throats connecting the pore bodies in a 10×20 region (termed trapped region hereafter and is delimited using a red box in Fig. A1) to the pore bodies outside of this region are selected so that their threshold pressure is larger than 1.4 bar. Simulation is then performed with an applied pressure difference of 1.4 bar across the entire pore network. The solvent liquids in the trapped region would be permanently confined in this region without evaporation effects. Figure A1. The evolution of the liquid saturation distribution inside the “trapped” region (marked using a red box), where a liquid cluster would be permanently trapped in absence of the evaporation effects. The evolution of the liquid saturation distribution inside the trapped region is shown in Fig. A1. Due to evaporation-induced emptying of pore throats surrounding the trapped region, purge gas breaks into the liquid cluster in this region at t ~ 95.0 ms and displaces the liquid solvent downstream (see the snapshot at point A in Fig. A1). The gas flow then fingers through the liquid cluster in the trapped region, which in turn generates new gas pathways (see the snapshot at point B). During this process, liquids are drained out of the trapped region. Meanwhile, the liquid drainage and gas flow fingering change the distribution of the remaining liquids in and around the trapped region (see the snapshot at point C), which greatly affects the evaporation of liquids in the partially-saturated pore bodies in the trapped region. Because of the coupled drainage and evaporation effects, the liquids in the trapped region are dried out eventually. Acknowledgements: The allocation of computing time by the ARC at Virginia Tech is acknowledged. Declarations of interest: none. References 1. Panagiotou, N.; Stubos, A.; Bamopoulos, G.; Maroulis, Z., Drying Kinetics of a Multicomponent Mixture of Organic Solvents. Drying technology 1999, 17, 2107-2122. 20 2. Peishi, C.; Pei, D. C., A Mathematical Model of Drying Processes. Int. J. Heat Mass Tran. 1989, 32, 297-310. 3. Bruin, S.; Luyben, K. C. A. In Drying of Food Materials: A Review of Recent Developments, 1st. Int. Symp. Drying, Montreal, 1978, 1978. 4. Le Gallo, Y.; Le Romancer, J.; Bourbiaux, B.; Fernandes, G. In Mass Transfer in Fractured Reservoirs During Gas Injection: Experimental and Numerical Modeling, SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers: 1997. 5. Sharma, G.; Verma, R.; Pathare, P., Thin-Layer Infrared Radiation Drying of Onion Slices. Journal of Food Engineering 2005, 67, 361-366. 6. Mujumdar, A. S., Handbook of Industrial Drying; CRC press, 2014. 7. Defraeye, T., Advanced Computational Modelling for Drying Processes–a Review. Applied Energy 2014, 131, 323-344. 8. Waananen, K.; Litchfield, J.; Okos, M., Classification of Drying Models for Porous Solids. Drying technology 1993, 11, 1-40. 9. Luikov, A. V., Heat and Mass Transfer in Capillary-Porous Bodies. In Advances in Heat Transfer, Elsevier: 1964; Vol. 1, pp 123-184. 10. Hirt, C. W.; Nichols, B. D., Volume of Fluid (Vof) Method for the Dynamics of Free Boundaries. J. Comput. Phys. 1981, 39, 201-225. 11. Prodanović, M.; Bryant, S. L., A Level Set Method for Determining Critical Curvatures for Drainage and Imbibition. Journal of colloid and interface science 2006, 304, 442-458. 12. Aidun, C. K.; Clausen, J. R., Lattice-Boltzmann Method for Complex Flows. Annual review of fluid mechanics 2010, 42, 439-472. 13. Philip, J.; De Vries, D., Moisture Movement in Porous Materials under Temperature Gradients. Eos, Transactions American Geophysical Union 1957, 38, 222-232. 14. Yiotis, A. G.; Stubos, A. K.; Boudouvis, A. G.; Tsimpanogiannis, I. N.; Yortsos, Y. C., Pore-Network Modeling of Isothermal Drying in Porous Media. Transport in Porous Media 2005, 58, 63-86. 15. Yiotis, A. G.; Boudouvis, A. G.; Stubos, A. K.; Tsimpanogiannis, I. N.; Yortsos, Y. C., Effect of Liquid Films on the Drying of Porous Media. Aiche Journal 2004, 50, 2721-2737. 16. Prat, M., On the Influence of Pore Shape, Contact Angle and Film Flows on Drying of Capillary Porous Media. International Journal of Heat and Mass Transfer 2007, 50, 1455-1468. 17. Plourde, F.; Prat, M., Pore Network Simulations of Drying of Capillary Porous Media. Influence of Thermal Gradients. International Journal of Heat and Mass Transfer 2003, 46, 1293-1307. 18. Wu, R.; Zhao, C.; Tsotsas, E.; Kharaghani, A., Convective Drying in Thin Hydrophobic Porous Media. Int. J. Heat Mass Tran. 2017, 112, 630-642. 19. Wu, R.; Kharaghani, A.; Tsotsas, E., Two-Phase Flow with Capillary Valve Effect in Porous Media. Chemical Engineering Science 2016, 139, 241-248. 20. Wilkinson, D.; Willemsen, J. F., Invasion Percolation: A New Form of Percolation Theory. Journal of Physics A: Mathematical and General 1983, 16, 3365. 21. Lenormand, R.; Zarcone, C.; Sarr, A., Mechanisms of the Displacement of One Fluid by Another in a Network of Capillary Ducts. J. Fluid Mech. 1983, 135, 337-353. 22. Yoon, R.-H. Methods of Enhancing Fine Particle Dewatering. Patent No. 6,855,260, 2005. 23. Xiumin, J.; Chuguang, Z.; Che, Y.; Dechang, L.; Jianrong, Q.; Jubin, L., Physical Structure and Combustion Properties of Super Fine Pulverized Coal Particle. Fuel 2002, 81, 793-797. 21 24. Joekar-Niasar, V.; Hassanizadeh, S. M.; Dahle, H., Non-Equilibrium Effects in Capillarity and Interfacial Area in Two-Phase Flow: Dynamic Pore-Network Modelling. J. Fluid Mech. 2010, 655, 3871. 25. Joekar-Niasar, V.; Hassanizadeh, S., Analysis of Fundamentals of Two-Phase Flow in Porous Media Using Dynamic Pore-Network Models: A Review. Critical reviews in environmental science and technology 2012, 42, 1895-1976. 26. Prat, M., Percolation Model of Drying under Isothermal Conditions in Porous Media. International Journal of Multiphase Flow 1993, 19, 691-704. 27. Gielen, T.; Hassanizadeh, S.; Leijnse, A.; Nordhaug, H., Dynamic Effects in Multiphase Flow: A PoreScale Network Approach. In Upscaling Multiphase Flow in Porous Media, Springer: 2005; pp 217236. 28. Thompson, K. E., Pore‐Scale Modeling of Fluid Transport in Disordered Fibrous Materials. Aiche. J. 2002, 48, 1369-1389. 29. Khayrat, K.; Jenny, P., A Multi-Scale Network Method for Two-Phase Flow in Porous Media. J. Comput. Phys. 2017, 342, 194-210. 30. Krishna, R.; Wesselingh, J., The Maxwell-Stefan Approach to Mass Transfer. Chemical Engineering Science 1997, 52, 861-911. 31. Mason, E. A.; Malinauskas, A., Gas Transport in Porous Media: The Dusty-Gas Model; Elsevier Science Ltd, 1983; Vol. 17. 32. Brown, G. P.; DiNardo, A.; Cheng, G. K.; Sherwood, T. K., The Flow of Gases in Pipes at Low Pressures. Journal of Applied Physics 1946, 17, 802-813. 33. Vignes, A., Diffusion in Binary Solutions. Variation of Diffusion Coefficient with Composition. Industrial & Engineering Chemistry Fundamentals 1966, 5, 189-199. 34. Buddenberg, J.; Wilke, C., Calculation of Gas Mixture Viscosities. Industrial & Engineering Chemistry 1949, 41, 1345-1347. 35. Young, L. C.; Finlayson, B. A., Mathematical Models of the Monolith Catalytic Converter: Part Ii. Application to Automobile Exhaust. Aiche. J. 1976, 22, 343-353. 36. Davis, T. A., Direct Methods for Sparse Linear Systems; Siam, 2006; Vol. 2. 37. Rothe, E. W.; Bernstein, R. B., Total Collision Cross Sections for the Interaction of Atomic Beams of Alkali Metals with Gases. J. Chem. Phys. 1959, 31, 1619-1627. 38. Liu, H.; Valocchi, A. J.; Kang, Q.; Werth, C., Pore-Scale Simulations of Gas Displacing Liquid in a Homogeneous Pore Network Using the Lattice Boltzmann Method. Transport. Porous. Med. 2013, 99, 555-580. 39. Homsy, G. M., Viscous Fingering in Porous Media. Annual review of fluid mechanics 1987, 19, 271311. 22