A numerical model for the transport and drying of solutions in thin porous media

We have developed a comprehensive numerical model for the transport and drying of solutions in thin porous media that consist of permeable fibers such as paper. We explicitly account for the gas-phase transport dynamics. Moreover, we introduce an empirical relation for the concentration-and molecular-weight dependence of the pore-fiber transport rate of the solutes. These two key elements enable us for the first time to realistically model two important phenomena relevant to inkjet printing technology. The first is the equivalent of the coffee-stain effect for dilute solutions in porous media. The second is the formation of solute rings for concentrated aqueous mixtures of compounds with a molecular weight significantly above that of water. Whereas the first is governed by spatially non-uniform solvent evaporation, the second case is dominated by solvent-mediated pore-fiber transport. We achieved a good qualitative agreement with the available experimental data.


Introduction
The formation of a coffee stain at the perimeter of an evaporating solution droplet resting with pinned contact lines on an impermeable surface is a well-known phenomenon [1].Dou & Derby, Huang et al. and Li et al. studied coffee-stain formation on top of porous media [2][3][4], for which no solute transport inside the porous media occurs.Coffee stains are usually not observed after deposition of dilute are typically water as an environmentally safe solvent, colorants, surfactants as well as so-called co-solvents [16].The latter are polar liquids of low volatility such as glycerol, with a concentration up to 50% [16].Their primary purpose is the prevention of inkjet nozzle clogging [16][17][18], among others [18][19][20].Wijburg et al. [21] investigated the transport of non-dilute, aqueous co-solvent solutions in paper substrates and observed the occurrence of solute rings after drying.
Our main interest is to understand the redistribution dynamics of solutes contained in water-based inks after printing, and specifically the formation of coffee-stains and solute rings inside thin porous media.For this purpose, we studied the transport processes after droplet deposition during the drying phase of paper sheets that have an initial maximum moisture content  comparable to the maximum holding capacity of paper, which is approximately 1.0 kg water/1 kg of dry paper [22].The redistribution of solid-phase solutes such as typical dyes is essentially arrested when the moisture content falls below a threshold of  crit ≈ 0.07.Below  crit there is essentially only adsorbed water left in the paper sheet, but the mobile liquid-phase water has evaporated.Therefore, we consider only liquid-phase ink transport in our model, i.e. a combination of Darcy flow and liquidphase diffusion and dispersion.This approach is obviously also suitable for solutes that remain in a liquid state after water evaporation such as typical co-solvents.Their redistribution is not arrested when most of the water has evaporated, but it markedly slows down due to the sensitive dependence of viscosity on solvent concentration and of the effective permeability on liquid content.
Paper is a thin porous medium consisting mainly of permeable cellulose fibers [23].Usually the presence of permeable fibers is not explicitly considered in the literature on paper drying [24][25][26][27][28][29].This is not necessary, because the timescale for pore-fiber equilibration, which is on order of 1 s for liquid-phase water in unsized paper, is typically much shorter than the drying time.In contrast, in the final stage of paper drying below  crit , this equilibration time is much longer (e.g.300 s in [30]).Therefore, Bandyopadhyay et al. presented a model, assuming water vapor diffusion is the dominant transport mechanism, that does take the presence of permeable fibers explicitly into account [30].Zapata et al. improved on this model by also considering heat transfer [31].They assumed the pore-fiber transport rate to be an empirical, linear relation with a constant time-scale, which implies that its validity is restricted to modest variations of the moisture content.
In our case, the initial moisture distribution is not spatially uniform, but highly heterogeneous, due to the localized deposition of an ink droplet on a paper sheet.Venditti et al. presented a model for the transport and evaporation of dilute solutions of surface-active compounds in paper.They took the presence of permeable fibers into account, although the pore-fiber equilibration time is short for dilute solutions [32].This is an essential element, however, because the fibers can take up both solvent and solute that become essentially immobilized temporarily, which greatly affects the transient moisture profile and the solute distribution upon droplet imbibition [32].
Venditti et al. approximated the evaporative solvent flux  ev by a global mass transfer relation where,  H 2 O ( = 0) and  H 2 O ( → ∞) are the water vapor concentrations at the paper surface  = 0 and in the ambient atmosphere, respectively, and  ∞ is a mass transfer coefficient.Expressions of this type are commonly used in the drying literature [24,25,[27][28][29]31,[33][34][35].Eq. ( 1) can properly account for the overall mass loss, but cannot reproduce its position dependence near boundaries of wet zones.Therefore, we adopt the diffusion-limited evaporation model used e.g. by Murali et al. [36], where  H 2 O is the diffusion coefficient of water vapor in air.We solve for flow, heat transfer and water vapor Viscosity  TEG and surface tension  TEG as a function of composition for aqueous TEG solutions [40].
transport in the gas phase, coupled to the flow, heat and mass transfer inside the paper sheet.This yields the correct position dependence of the evaporative flux, even at abrupt wet zone boundaries, which is a key ingredient to study coffee-stain and solute-ring formation.
To account for the micro-structure of paper, we applied a dualmedium description [32], which allows for transport of liquid and solutes between the inter-fiber and intra-fiber porespaces.Typical cosolvents have a molecular weight (MW) between 50 and 300.If their MW is below the cut-off value for penetration of the cellulose fiber walls [37][38][39], they can migrate into the cellulose fibers.Wijburg et al. found that this pore-fiber migration rate is highly sensitive to the cosolvent MW and the local co-solvent concentration [21].They ascribed the latter effect to water acting as a plasticizer of the cellulose fiber walls.We extended the pore-fiber exchange model introduced in [32] to account for both this plasticization effect and the dependence on solute MW.We performed a systematic parameter study and obtained good qualitative agreement between our simulations and the experiments regarding coffee-stain formation by Nilghaz et al. [9] as well as the co-solvent ring formation experiments in [21].

Material properties
The material properties of mixtures and solutions are generally concentration dependent.For the case of aqueous co-solvent solutions, the parameters that are most sensitive to composition variations are surface tension, viscosity and saturated vapor pressure.For most of our simulations, we use the properties of tetra(ethylene glycol) (TEG)water mixtures as a representative co-solvent solution.Fig. 1 shows the viscosity  TEG and surface tension  TEG of TEG solutions, which vary by factors of approximately 50 and 0.6, respectively, over the composition range.Analytical representations of  TEG and  TEG are provided in the Supplementary Information, along with material properties of other common co-solvents.Many properties are only weakly composition dependent, such that we substitute values for pure substances for convenience.These approximations are also listed in the Supplementary Information.

Water activity
Both paper and co-solvents are hygroscopic materials that absorb appreciable quantities of water depending on the relative humidity  ℎ of the surrounding gas phase [41,42].For binary mixtures of water and typical co-solvents, the Van Laar relation [43,44] ln can represent the equilibrium mixture composition as a function of  ℎ quite well, provided the coefficients  12 and  21 are treated as adjustable parameters.Here,   is the solvent activity coefficient [45] and  cs and   = 1 −  cs are the mole fractions of co-solvent and water, respectively.Fig. 2 shows an example for TEG with fitted values  12 = −2.4 and  21 = −0.6.Analogous fits for selected other co-solvents are given in the Supplementary Information.
In equilibrium, the co-solvent weight fraction where   and  cs represent the mixture masses of water and co-solvent, respectively, can be translated into a moisture content (with units kg of water/kg of co-solvent) which is plotted as the black solid line in Fig. 2(b).It follows from Eq. ( 5) that For paper, the Guggenheim-Anderson-De Boer (GAB) isotherm can represent the moisture sorption behavior as a function of water activity well.Here,  0 ,  and  are adjustable parameters,  H 2 O is the water vapor mass density in the gas phase and  sat its saturated value.A typical example of a GAB isotherm is shown as the blue solid line in Fig. 2

(b).
To the best of our knowledge, the moisture isotherms of paper sheets containing various quantities of co-solvents per unit area are not available in the literature.Therefore, we make two assumption to arrive at a relation between water activity, solution mass and composition.Firstly, we assume that as far as their moisture sorption behavior is concerned, paper and co-solvent do not interact, i.e. the moisture sorption capacity of paper is not affected by the presence of and contact with co-solvents and vice versa.Secondly, we assume that the chemical potential of water is uniform, i.e. it has the same value in both the paper phase and the solution phase.
When solving the dynamic equations for liquid transport in paper, the quantity of co-solvent and water per unit area,  tot,cs and  tot,w , are known.For solving the gas-phase dynamics, the water activity at the paper-gas interface is required.Eq. ( 7) can be inverted [46,47] to yield 2 + ( The condition of local thermodynamic equilibrium then can be formulated as (10) and correspondingly the total local areal density of moisture is assigned to the paper phase and the co-solvent phase according to Here,  paper is the grammage or areal mass density of the porous medium.Thus, Eqs.(10), (11) represent two equations for the two unknowns  paper and  cs , which we solved numerically and from which   as a function of  tot,cs and  tot,w can be determined using Eq.(10).

Liquid transport in paper
We adopted and extended the model of [32], rewritten for an axisymmetric system.All quantities are assumed uniform in the thickness direction of the thin porous medium.This implicitly assumes that the timescale of transport in the thickness direction is much shorter than for the in-plane direction, which is consistent with the experimental systems we will compare to.The mass conservation equation for the liquid in the inter-fiber pores of the paper substrate is given by where  denotes the liquid content in the inter-fiber pores (units kg/m 2 ),   the substrate thickness,  liq the mass density of the liquid mixture,  the radial Darcy velocity,  , the evaporative flux from inter-fiber pores and   the pore-fiber transport rate, respectively.The radial Darcy velocity  is given by For the effective permeability  of the porous medium, we adopt the Van Genuchten relation as detailed in [22,32].Similarly, we adopt a Van Genuchten relation for the capillary pressure , which is determined by  and the solute concentration in the pores   =   () Here,  max denotes the maximum holding capacity of the inter-fiber pores and   and   are constants.The prefactor   () is a function of surface tension   and is thus composition dependent [32,48,49]   () ≡   (0)   () cos ()   ( = 0) cos ( = 0) , (15) where   is the liquid-air interfacial tension of the mixture and  is the contact angle on the pore walls.
We assume the permeability of the intra-fiber pores to be so small that the liquid can be considered as effectively immobile.Hence, the liquid content in the fibers   only changes due to the evaporative flux  , and pore-fiber transport We assume the capillary pressure in the intra-fiber pores to be higher than in the inter-fiber pores: Here,  f ,max denotes the maximum holding capacity of the intra-fiber pores and   is a constant.For the commercial, inkjet-optimized printing paper used in the experiments in [21], it holds to good approximation that  max =  f,max and that () =   () = 0. Fig. 3(a) illustrates the capillary pressures in the pores () and in the fibers   (  ) for different values of   .In thermodynamic equilibrium, the capillary pressures of the liquid in the inter-fiber and intra-fiber pores are equal.Consequently, the parameter   determines the pore-fiber distributions of the liquid in the equilibrium state.For   = 1,  and   are equal in equilibrium.This case corresponds to the blue, solid, diagonal line in Fig. 3(b), where the non-dimensional equilibrium liquid content in pores ∕ max (solid lines) and fibers   ∕ f,max (dashed lines) are plotted as a function of the total liquid content  tot ≡  +   for different values of   .For   > 1, the equilibrium content corresponding to a specific value of the capillary pressure will be higher in the fibers than in the pores.For   → ∞ and for increasing total liquid contents  tot , the pores will remain essentially empty, while the fibers take up essentially all the liquid until the fiber holding capacity is reached.Only beyond this point, i.e. for  tot >  f,max , will the pores start to fill, which is illustrated by the green solid and dashed lines in Fig. 3(b).

Pore-fiber transport
In [32], we considered dilute solutions and assumed the pore-fiber transport rate   to be proportional to the pressure difference   (  )− () between pores and fibers [50].The MW of co-solvents is typically much higher than that of water (see Tab. 1 in the Supplementary Information), which tends to increase viscous friction and reduce   .Moreover, in [21] we concluded that the presence of a sufficient quantity of water swells and plasticizes the fiber walls and thereby facilitates the pore-fiber transport of bigger molecules.To account for these effects, we propose to modify   as follows Here,   ( ) is an MW-dependent prefactor.The solution viscosity ( , ) is strongly concentration dependent, see Fig. 1(a).The function (  ) represents the effective permeability of the cellulose fiber walls.We assume it to be given by where  = 0.02.The function  () represents the influence of the co-solvent concentration in the pores.We assume it to be given by: where   = 0.005.The parameter  crit denotes the threshold above which the co-solvent concentration starts to affect the pore-fiber transport rate.The parameters   and   determine the abruptness of the decay of   beyond  crit and the baseline rate for  ≫  crit , respectively.Fig. 4 illustrates   for  crit = 470 kg∕m 3 and three values of   .

Solute and heat transport in paper
We adopt the model and notation of [32] and solve Eqs.(16,(23)(24)(25)) and ( 27) in [32], rewritten for an axisymmetric system, for the time evolution of the solute concentration in the inter-fiber pores  and the intra-fiber pores   (units kg of solute per m 3 of solution), the adsorbed solute concentration  ad as well as the paper temperature  paper .The terms   ,   and   are set to zero.In the following we only point out the modifications and extensions with respect to the model detailed in [22,32].
For axisymmetric flow in the radial direction, only longitudinal dispersion occurs.We use the following expression for the dispersion coefficient  ef f [51][52][53]  where the local Peclet number  is given by where  pore = 10 μm is a typical longitudinal lengthscale (cellulose fiber width) and   is the molecular diffusion coefficient.Eq. ( 21) has been obtained for random sphere packings.To the best of our knowledge, experimental or computational studies for fibrous, paper-like materials are not available.
For the transport of co-solvent solutions, adsorption effects are usually weak, because the adsorption capacity is much smaller than the deposited co-solvent content.In contrast, for the transport of dilute solutions, adsorption is a strong effect, contributing to chromatographic retention [5].

Gas phase transport
We adopt the model of [36], rewritten for an axisymmetric system, and solve for the velocity ⃗  = (  ,   ) and the temperature  gas of the gas phase, as well as the water vapor mass density  H 2 O .

Coupling the paper and gas domains
In contrast to [22,32], the heat exchange with the gas phase via thermal conduction and evaporative cooling is now accurately accounted for by the expression where   is the enthalpy of evaporation of water and   is the evaporative mass flux.Specifically, the position dependence of the evaporative flux is accurately accounted for by the expression where  H 2 O is the diffusion coefficient of water vapor in air.The evaporative flux is split into a contribution from the inter-fiber pores  , and a separate one from the intra-fiber pores  , Here,   = 0.2 is a constant that accounts for the increased mass transfer resistance for evaporation from the fibers, as the water first has to go through the fiber walls.In Eqs. ( 23), (24) we only consider the gas phase above the paper for simplicity.
The coupling between the paper-and the gas domain is done by assuming continuity of temperature Moreover, the water vapor density  H2O at the paper surface  = 0 is determined by the effective sorption isotherm of the porous medium as described in Section 2.1

Boundary conditions (BCs)
All simulations are conducted assuming axisymmetry, with a gas domain of radius  = 0.1 m and height  = 0.2 m.The thermodynamic variables on the upper and right boundaries of the gas phase domain are assumed to be determined by the ambient conditions (temperature  amb = 298 K, relative humidity  ℎ , atmospheric pressure  amb = 1 bar) Moreover, we assume the radial and vertical velocity components (  ,   ) to be zero there as well as on the lower boundary  = 0 The BCs on the left boundary  = 0 represent axisymmetry.An overview of all BCs is provided in Fig. 5 for the case where the radial extension of the porous substrate   is smaller than .The BCs for the dynamics inside the porous medium at its perimeter are no-flux conditions, i.e. all radial derivatives vanish at  = 0 and  =   .
We have also considered systems with   =  or where a gas phase is present also below the paper substrate.An adapted version of Fig. 5 is presented in the Supplementary Information.

Initial conditions (ICs)
We do not model the deposition of the liquid.Rather, we assume that the imbibition process is already complete at  = 0 and that an Here,  0 is the initial solute concentration,  0 the initial deposited ink content and f hs is a smooth Heaviside function with two continuous derivatives The parameter  denotes the transition width, i.e. a smaller value of  corresponds to a more abrupt transition.We assume that no solute is present in the fibers initially, i.e.   ( = 0) = 0, but that the fibers contain the same quantity of water   (amb) as dry paper in equilibrium with the ambient atmosphere.Similarly, the same quantity of water (amb) is contained in the pores at  >  0 +  as present in dry paper in equilibrium with the ambient atmosphere.Both   (amb) and (amb) are determined by the effective sorption isotherm of the porous medium and the condition of equal capillary pressure.
We assume there is no initially adsorbed solute We initialize the gas phase with zero velocity and constant temperature If we would assume the gas phase to have an initially uniform water vapor concentration of  H 2 O (amb) ≡  ℎ  sat ( amb ), then the initial evaporation rate would be infinite above the wet zone, which is numerically problematic.Therefore, we assume that  H 2 O is initialized with an exponentially decaying function above the wet zone The amplitude () is equal the difference between  H 2 O (amb) and Eq. ( 28); the decay length  0 ≈ 1 mm.The influence of this IC disappears typically within a few seconds, which is short compared to the drying time, which is typically 5-10 min.

Numerical results
We performed a systematic parameter study using the material parameters of TEG.Unless mentioned otherwise the default values of important parameters are those listed in Table 1.A complete list of all symbols and parameter values is provided at the end of the Supplementary Information.

Model validation
For the purpose of validation we considered the drying dynamics of a papersheet after deposition of a droplet of pure water, for which quantitative experimental data were available in [21].Details are provided in the Supplementary Information.The gas density  gas is affected by both temperature and humidity variations, which induce buoyancy-driven convection.Where water evaporates from paper, the local humidity in the adjacent gas increases, but the temperature drops due to evaporative cooling.As  gas decreases with increasing humidity but increases with decreasing temperature, these two mechanisms are antagonistic.Temperature variations are on order of a few K corresponding to a relative density change  gas ∕ gas ≈ 1%.Variations in relative humidity are on order of 60%, corresponding to  gas ∕ gas ≈ 2%.Therefore, humidity variations tend to govern the large scale convection pattern.

Gas phase dynamics
In Fig. 6(b), a plume with upwards flow is observed above the wet zone.Its maximum velocity is about |⃗ | max ≈ 1 cm/s.This corresponds to a Reynolds number   =  gas  max ∕ gas ≈ 60, which implies that inertial effects are important, but that the flow is laminar.Here  = 10 cm is the width of the computational domain.The value of the Peclet number is   =  max  0 ∕ H2O = 2, which indicates that the water vapor transport is essentially diffusion dominated, but that convection effects cannot be neglected.The same conclusion holds for heat transfer in the gas phase, as the Lewis number of air  =  air ∕ H2O is close to 1. Here,  air is the thermal diffusivity of air.radii of the wet zones.Fig. 7(b) shows the radius of the wet zone  wet zone and the width of the ring stain  ring (defined as the full width at half maximum) as a function of  .The blue solid line corresponds to a fit function

Co-solvent transport
with fitparameters ,  and  0 .The exponent 1∕2 follows from mass conservation for drop volumes sufficiently large to completely penetrate the substrate thickness   .The green dashed line in Fig. 7(b) represents a constant value of  ring .Both curves are in good qualitative agreement with Fig. 5(b) in [21].However, in the experiments in [21], the amplitude of the ring stain decreases with decreasing volume, whereas in Fig. 7(a) it increases.In the three curves in Fig. 7(a), the value of  is identical.It is conceivable that in the experiments  was not independent of  .Fig. 8(a) shows the co-solvent distribution in the inter-fiber pores for different values of the transition width  at  = 160 s.The amplitude of the ring stain decreases with increasing .The reason is that the local water evaporation rate is higher for sharper transitions.As a result, the concentration  exceeds the threshold  crit sooner and thus more co-solvent is trapped in the pores.Since a sharper initial distribution causes a stronger Darcy flow, the radius of the wet zone increases slightly more for smaller value of .
Fig. 8(b) shows co-solvent content in the pores for different values of   .Larger values of   reduce the quasi-steady co-solvent content in the inter-fiber pores.There are two effects at play.Firstly,   increases the pressure differential in Eq. ( 18), i.e. the driving force for pore-fiber transport.A larger value of   thus reduces the timescale for porefiber equilibration and decreases the co-solvent content in the pores.A second factor is the higher equilibrium co-solvent content in the interfiber pores for lower values of   , as illustrated by Fig. 3.This will lead to considerable co-solvent contents in the inter-fiber pores -even if pore-fiber equilibrium is reached.
As shown in Fig. 8(b), the radius of the wet zone increases for smaller value of   .The reason is that smaller value of   implies slower pore-fiber transport.Consequently, Darcy flow induced spreading occurs at a higher rate, primarily because the permeability  is a nonlinear function of  [22].
Fig. 8(c) shows the co-solvent content in the inter-fiber pores after a quasi-steady state is obtained.The maximum amplitude of the cosolvent content in the ring-stain increases for smaller   .This is because a faster decay of the function  () reduces the pore-fiber transport rate sooner, and thus leaving more co-solvent trapped in the inter-fiber pores.
Fig. 8(d) shows co-solvent distribution for different   at  = 160 s.For smaller   , liquid transport between pores and fibers is slower, whereas evaporation is unchanged.As a result, more co-solvent will be trapped in the pores, especially near the perimeter of the wet zone, where evaporation is highest.The larger amount of liquid in the fibers, the smaller will be the wet-zone diameter, because liquid is assumed immobile in the fibers according to Eq. ( 16).For   = 7 ⋅ 10 −9 and 3.5 ⋅ 10 −9 kg/m 2 , a pronounce ring stain is observed.Its amplitude and width are larger for the smaller value of   .This is in good qualitative agreement with the experimental data for TEG and PEG-6 in Fig. 3(g) of [21].Fig. 8(e) shows the co-solvent distribution in the inter-fiber pores for different initial concentrations  0 at  = 160 s.For larger values of  0 , more co-solvent is retained in the inter-fiber pores.The reason is that the threshold  crit is exceeded from the outset for  0 ≥ 50 wt%.
In Fig. 8(f) we plotted the total co-solvent content in pores and fibers  cs +  cs,f for the default parameter values in Table 1.It is clearly observed that no pronounced peak is visible at the perimeter of the wet zone, in contrast to the peaks in the inter-fiber pore content observed in e.g.Figs.7(a) or 8(a,b).Rather, the overall distribution resembles that of an initially uniform distribution that has undergone some radial spreading in a regular porous medium with impermeable solid phase.Only a small shoulder is present at  ≈ .The shape qualitatively agrees with the experimental data for the capacitance increment in Fig. 6(b) of [21], which is sensitive to the total liquid content.

Evaporative flux near wet zone boundaries
The coffee stain effect for sessile droplets on impermeable surfaces is enhanced by the strong non-uniformity of the evaporative flux, which exhibits a singularity of the form for the case of small contact angles and slow, diffusion-limited evaporation [1].Here,  drop is the footprint radius of the droplet, which is assumed to have a spherical-cap shape.The singularity arises from the abrupt change in BC for solvent vapor transport from a constant vapor pressure at the liquid-air interface to a no-flux condition beyond the droplet footprint.Although unphysical, the singularity is integrable, i.e. the overall mass flux remains finite.
For slow, diffusion-limited evaporation, the time-derivative in Eq. ( 12) of [36] can be neglected [54].Moreover, near the surface at  = 0, the air velocity is small owing to the no-slip BC.Therefore, the convective term in Eq. ( 12) of [36] can be neglected as well, leading to the Laplace equation After droplet deposition on a porous medium, there is no abrupt change in the type of BC, which for volatile droplets on impermeable substrates changes from Dirichlet-to Neumann type.However, there can still occur an abrupt change in the magnitude of the vapor pressure according to Eq. ( 8), which is a Dirichlet BC for Eq. ( 41).In the next paragraph, we show that this also leads to a singularity in the evaporative flux.Using an integral transform approach [55], the following analytical solution of the Laplace equation  = 0 for a semi-infinite 2D domain with piece-wise constant but discontinuous Dirichlet BCs can be derived Thus, the evaporative flux  ev () is proportional to as a function of || with Eq. (44).The latter exhibits an unphysical, nonintegrable singularity due to the boundary discontinuity.However, this does not pose an essential problem in practice, because in a homogeneous porous medium, a discontinuous liquid distribution would induce a singular Darcy flow, thus removing the discontinuity in an instant.In contrast, the (integrable) singular evaporative flux near the contact line of a sessile droplet persists in time, because the discontinuity of the liquid distribution remains present until the contact line depins or the evaporation is complete.(42).The symbols represent a numerical finite-element solution, the black solid line corresponds to Eq. ( 44).The dotted vertical lines indicate the mesh size in the vicinity of the discontinuity at  = 0 and the halfwidth of the simulation domain, respectively.The insets illustrate the BCs and the non-uniform finite-element mesh around the origin.

Coffee stain formation in thin porous media
After deposition of a droplet on a thin porous medium, the nonuniform liquid distribution sets up a non-uniform evaporative flux and induces a Darcy flux directed radially outward [see e.g.Fig. 5(b) in the Supplementary Information].Therefore, in principle, a coffee-stain-like effect can also occur for solutions in porous media.However, for dilute solutions, chromatographic retention [5] often depletes the solute near the perimeter of the wet zone.Moreover, the Darcy flow is usually decelerating, as the gradients in the liquid content become shallower in time [see e.g.Fig. 5(b) in the Supplementary Information], in contrast to the flow in a sessile volatile droplet, which accelerates in time [1].
Exceptions of practical relevance exist for wet zones that extend to the perimeter of a paper substrate or a hydrophobic barrier.In these cases, the liquid distribution exhibits a persistent discontinuity at the edges of the porous sheet or the wettable region.The flux singularity is then expected to be integrable with exponent −1∕2, due to the no-flux condition relevant beyond the perimeter of the paper sheet or on top of the hydrophobic barrier [9][10][11][12][13][14][15].Fig. 10(a) shows numerical simulations of the solute distribution  cs in the inter-fiber pores for circular paper substrates with finite radius   that were initially homogeneously wet.The ICs used deviate slightly from those in Section 3.7 and Table 1 in that A pronounced peak in the solute content is observed at the substrate edge at  =   = 5 mm for sufficiently small values of   .Fig. 10(b) presents a zoom of (a) near  =   .Fig. 10(c) shows the width of the solute peak  as a function of   .The width is defined as the halfwidth-at-half-maximum (HWHM) between the minimum of  cs at  = 0 and its maximum at  =   , as indicated by the red arrow in Fig. 10(a).The red solid line in Fig. 10(c) is a power law fit  ∼  0.55  .Nilghaz et al. presented experiments of the evaporation of dye solutions in paper, and observed pronounced coffee-stains at hydrophobic barriers and edges of finite-area paper substrates (see Figs. 7 and 9(b) in [9]), whereas no coffee stains were observed for droplets deposited on large-area substrates.Our simulations in Fig. 10 qualitatively reproduce their results.Fig. 10(c) suggests that an estimate for the solute diffusion coefficient can be obtained from the width of the coffee stain.

Difference between coffee stains & solute rings
The distinguishing feature between the co-solvent rings observed in Ref. [21] and e.g.Fig. 8(d) and the coffee-stains in Ref. [9] and Fig. 10, is that the latter involve sustained solute transport towards the edge.This is illustrated in Fig. 10(d), where the normalized Darcy velocity is plotted at fixed radial positions  * .The solid line represents the simulation shown as the red curve in Fig. 7 at  * = ∕ 0 = 1.The dashed line corresponds to the curves in Fig. 10(a) at  * = ∕ 0 = 0.9.While the velocity decays to 10% of its initial value within approximately 40 s for the red curve, it remains essentially constant for the blue one, during the time interval shown.Moreover, while the co-solvent content in the ring-stains observed in e.g.Fig. 7 is lower than the initial cosolvent content  0  0 in the inter-fiber pores, the solute content in the coffee-stains in Fig. 10(a,b) is higher than the initial value for   ≤ 10 −10 m 2 /s.

Pore-fiber transport -effect of solute concentration and molecular weight
The selected of  crit = 470 kg/m 3 is motivated by the experiments in [21], where it was observed the residual co-solvent content left in the pores of a paper substrate after deposition and evaporation of an aqueous TEG solution droplet noticeably increased beyond an initial concentration  0 of approximately 50 wt%.Moreover, it was reported in Ref. [56] that the persistent expansion strain after deposition and evaporation of various aqueous co-solvent solutions monotonically increased for low co-solvent concentrations up to a critical value between 40 and 60 wt%, beyond which it monotonically decreased.This critical value did not seem to depend strongly on the chemical composition nor the MW of the co-solvents considered in Refs.[21,56].The co-solvent MW therefore appears to primarily affect the values of   and possibly also of   and   , but not so much  crit .
Peppas and Reinhart studied solute transport through swollen polymer membranes [57].They proposed a linear correlation between the logarithm of the membrane permeability for the solute and the square of the solute molecular radius  2  .For solutes of low MW,   is expected to scale as ( ) 1∕3 .Yasuda et al. [58] and Sato & Kim [59] reported experimental data consistent with this prediction.
Tokita [60] studied the diffusion of probe molecules in cross-linked polyacrylamide gel and observed a good correlation between the gel diffusion coefficient of the solutes and ∕ 0 = exp(−) , where  ∼ ( ) 1∕3  −3∕4 tot ,  0 is the diffusion coefficient in pure solvent and  tot is the mass density of the gel (not including the solvent) in kg/m 3 .This result is similar to the relation proposed by Ogston et al. [61] and Cukier [62], who studied the diffusion of Brownian spheres in semidilute polymer solutions theoretically.It is consistent with the experimental results of Colton et al. [63] for probe diffusion in cellulosic membranes and of Langevin & Rondelez [64] for colloid transport through semidilute polymer solutions.In contrast, Brown and Johnsen [65], who studied the MW dependence of the effective diffusion coefficient for low MW solutes through polyacrylamide gels, found a correlation  ∼

𝑠
for an agarose gel and   ≲ 15 nm.The above experimental and theoretical results indicate that the pore-fiber transport rate is strongly dependent on MW.According to Eq. ( 18), the MW-dependent terms in the pore-fiber transport rate   are given by   ∕, i.e. both   and the solution viscosity  can contribute to the MW dependence of   .For non-dilute, aqueous solutions of poly(ethylene glycols), the scaling of  with MW is well represented by a powerlaw relation  ∼ ( )  with exponent 0.75 <  < 1, at least up to MW = 300 (see Fig. 3 in the Supplementary Information).
In Fig. 11(a) experimental data for the differential light transmittance ∕ 0 are shown, which are taken from Fig. 3(g) in [21].The three curves represent aqueous solutions of ethylene glycol (EG), TEG and hexa(ethylene glycol) (PEG6).The relation between ∕ 0 and  is approximately linear for sufficiently large , but non-linear for small  [21].In the corresponding numerical simulations shown in Fig. 11(b), we assumed a constant value of   = 2.4 ⋅ 10 −9 kg/m 2 for the different co-solvents EG, TEG and PEG6, i.e. that only the viscosity contributes to the MW dependence of   .The simulations resemble the experimental curves qualitatively well as for the ratio of the peak heights for TEG and PEG-6.Moreover, the width of the peaks increases with increasing MW in a comparable fashion in both Figs.11(a) and (b).In Fig. 11(c) we plotted the peak values of co-solvent pore content max( cs ) as a function of co-solvent MW.The solid line represents a linear relation, in excellent qualitative agreement with Fig. 3(h) in [21].

Rehydration
It was found in [21] that the deposition of pure water on regions where co-solvent is temporarily 'trapped' in the pores greatly facilitates pore-fiber transport.Fig. 12(a) shows an optical transmission image of a paper samples after deposition and drying of a 50 wt% TEG droplet.The gray area with increased transmission is due to residual co-solvent residing in the pores.Fig. 12(b) shows the same sample after deposition and evaporation of two small water droplets in the center.The black in the middle exhibits the same transmission as the dry paper in the periphery and arises a consequence of the completed, solvent-assisted pore-fiber transport of TEG [21].
We reproduce this process numerically by application of the following ICs: Here,  out = 0.2 max denotes the temporarily trapped co-solvent content in the pores and  in = 0.7 max the increased liquid content due to deposition of the rehydrating water distribution with radius   = 2 mm and abruptness parameter  = 0.5 mm.The initial co-solvent content in the pores  cs ( = 0) is assumed independent of .Thus, the co-solvent concentration  satisfies Eq. ( 51).Eq. ( 52) expresses that no co-solvent is initially present in the fibers.Figs.12(c,d) show the initial and the final co-solvent distribution in the pores, after the rehydrating water has evaporated.The similarity with Fig. 12(a,b) is striking.

Solutocapillary effects
Solutocapillary effects due to the concentration dependence of surface tension () are relatively weak in the co-solvent ring formation process, because they are essentially quenched by the strong increase of viscosity ().This is in stark contrast to dilute surfactant solutions [32], where () strongly influences the ink redistribution and () does not change as much.

Material parameters and closure relations
Models for the transport of liquids and solutions in unsaturated porous media based on the Richards equation are conceptually attractive, because they explicitly take the governing physical processes and corresponding state variables into account.Their main disadvantage is that they require a large number of material parameters and closure relations that are often not available from published experimental data.The fact that the detailed composition and properties of commercial paper grades are usually maintained confidential adds to the problem.Therefore, we resolved to using empirical relations, some of which are well-established in literature and some were newly introduced.Examples include the 1.Van Genuchten relations for the capillary pressures  and   in Eqs. ( 14), (17) and the effective permeability as a function of liquid content.2. Dispersion relation according to Eq. ( 21). 3. Langmuir solute adsorption dynamics according to Eq. (27) in [32].4. Water activity: we assume that Eqs. ( 3), (7), which are accurate for the two-component systems water/paper and water/cosolvent, also hold for the three-component system water/paper/ co-solvent.
To improve the accuracy of our model, the parameters relevant to items 1-4 in the above list need to be calibrated against dedicated experiments.The last two items would benefit from pore-resolved models and experiments [68][69][70][71][72] and dedicated single-fiber experiments [73][74][75][76].

Summary and conclusions
We have developed a numerical model for the transport and drying of solutions with volatile solvents in thin porous media.Our model explicitly includes the gas phase above the porous medium and therefore yields realistic evaporative flux distributions for the case of diffusionlimited evaporation.This is crucial for studying the coffee-stain effect in porous media, which is well reproduced by our model.In the case of sessile droplets on impermeable surfaces, the evaporative flux exhibits a singularity at the contact line with an exponent of −1∕2.In contrast, we found a flux singularity with exponent −1 for abrupt changes in liquid content in a homogeneous thin porous medium.In experiments, coffee-stains are primarily observed at the perimeter of media of finite dimensions or of moisture distributions confined by hydrophobic barriers.In these cases, the corresponding flux singularity exponent is −1∕2.
We also investigated the occurrence of solute rings deposition of concentrated aqueous solutions of co-solvents such as glycerol.Although their appearance resembles coffee-stains, their origin is different and associated with delayed pore-fiber transport in media with permeable fibers.We introduced an empirical relation for the concentrationand molecular-weight dependence of the pore-fiber transport rate of co-solvents.Our results indicate that the main contribution to the molecular-weight dependence stems from the solution viscosity.
Our model reproduces the available experimental data qualitatively well.It can therefore aid in the further development of inkjet printing technology, where transport and distribution of e.g.colorants needs to be tightly controlled to safeguard the color intensity and uniformity of printed patterns.Moreover, it can assist in the design of sensor devices based on paper microfluidics that e.g.capitalize on the concentration enhancement of solutes due to the coffee-stain effect.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig.1.Viscosity  TEG and surface tension  TEG as a function of composition for aqueous TEG solutions[40].

Fig. 2 .
Fig. 2. (a) Equilibrium co-solvent concentration of an aqueous solution of TEG as a function of relative humidity.The filled red circles and the open squares represent experimental data from Refs.[41,42], respectively.The solid line is a fit according to Eq. (3).(b) Moisture isotherms for TEG (black solid line) and paper (blue solid line, GAB parameters  0 = 0.04,  = 0.75 and  = 10).

Fig. 3 .
Fig. 3. (a) Capillary pressures in pores  (solid line) and fibers   (dashed lines) as functions of corresponding liquid contents  and   for different values of   .The vertical dotted line denotes the holding capacities  max =  f ,max .(b) Non-dimensional equilibrium liquid content in pores ∕ max (solid lines) and fibers   ∕ f ,max (dashed lines) as a function of the total liquid content for different values of   .

Fig. 4 .
Fig. 4. The function  (), governing the concentration dependence of the pore-fiber exchange rate   , for different values of   .

Fig. 5 .
Fig.5.Boundary conditions for the gas phase domain and the porous medium.The symbol   is short-hand for ∕.

Fig. 6 .
Fig. 6.Pseudo-color plots of (a) the normalized water vapor density  H 2 O ∕ sat ( amb ) and (b) the velocity distribution |⃗  gas | in the gas phase at  = 120 s after deposition of a droplet of pure water ( = 5 μl).The solid lines in (b) represent streamlines.

Fig. 6 (
Fig. 6(a) shows the normalized water vapor density  H 2 O ∕ sat ( amb ) above a paper substrate at  = 120 s after deposition of pure water, for an ambient relative humidity  ℎ = 29.7%.Fig. 6(b) shows the magnitude of gas velocity distribution and selected streamlines.The gas density  gas is affected by both temperature and humidity variations, which induce buoyancy-driven convection.Where water evaporates from paper, the local humidity in the adjacent gas increases, but the temperature drops due to evaporative cooling.As  gas decreases with increasing humidity but increases with decreasing temperature, these two mechanisms are antagonistic.Temperature variations are on order of a few K corresponding to a relative density change  gas ∕ gas ≈ 1%.Variations in relative humidity are on order of 60%, corresponding to  gas ∕ gas ≈ 2%.Therefore, humidity variations tend to govern the large scale convection pattern.In Fig.6(b), a plume with upwards flow is observed above the wet zone.Its maximum velocity is about |⃗ | max ≈ 1 cm/s.This corresponds to a Reynolds number   =  gas  max ∕ gas ≈ 60, which implies that inertial effects are important, but that the flow is laminar.Here  = 10 cm is the width of the computational domain.The value of the Peclet number is   =  max  0 ∕ H2O = 2, which indicates that the water vapor transport is essentially diffusion dominated, but that convection effects cannot be neglected.The same conclusion holds for heat transfer in the gas phase, as the Lewis number of air  =  air ∕ H2O is close to 1. Here,  air is the thermal diffusivity of air.

Fig. 7 (
Fig.7(a) shows the distributions of the co-solvent content in the inter-fiber pores  cs ≡ ∕ liq for different droplet volumes  , at  = 160 s.At this point in time, the co-solvent pore-fiber distributions are already quasi-steady.The vertical dashed lines denote the initial

Fig. 7 .
Fig. 7. (a) Co-solvent content in inter-fiber pores  cs for different deposited liquid volumes   = 160 s.(b) Radius of the wet zone  wet zone (open circles and blue solid line) and width of the ring stain  ring (green circles and green dashed line) as a function of  .

Fig. 8 .
Fig. 8. Co-solvent distribution in inter-fiber pores  cs for different values of (a) , (b)   , (c)   , (d)   and (e)  0 , all at  = 160 s.(f) Total co-solvent content  cs +  cs,f in pores and fibers for the default parameter values in Table1.

Fig. 9 .
Fig. 9. Boundary flux | | |   ( = 0) | | | as a function of || corresponding to the solution of the Laplace equation on a 2D halfspace with discontinuous Dirichlet boundary conditions according to Eqs.(42).The symbols represent a numerical finite-element solution, the black solid line corresponds to Eq. (44).The dotted vertical lines indicate the mesh size in the vicinity of the discontinuity at  = 0 and the halfwidth of the simulation domain, respectively.The insets illustrate the BCs and the non-uniform finite-element mesh around the origin.

Fig. 10 .
Fig. 10.(a) Solute distribution in the inter-fiber pores  cs for different values of   , for homogeneously wetted substrates with limited area.(b) Zoom of (a) near the edge of the substrate.(c) Halfwidth  of the solute peak near the substrate edge as a function of   .(d) Normalized Darcy velocity v( * , ) at fixed radial positions.The red solid line represents the red curve in Fig. 7 at  * = ∕ 0 = 1.The blue dashed line corresponds to the curves in (a) at  * = ∕ 0 = 0.9.

Fig. 11 .
Fig. 11.(a) Experimental data for differential light transmittance ∕ 0 after deposition and drying of aqueous EG, TEG and PEG6 solution droplets ( 0 = 40 wt%) in paper, taken from Fig. 3(g) in [21].(b) Numerical simulations of the residual co-solvent content in the pores after deposition and drying of aqueous EG, TEG and PEG6 solution droplets.Here, we have assumed   = 4.5,   = 2.4 ⋅ 10 −9 kg/m 2 ,  = 500 μm.(c) Peak values of co-solvent pore content as indicated by the horizontal arrows in (b) as a function of co-solvent MW.The solid line is a guide to the eye.

Fig. 12 .
Fig. 12. (a,b) Optical transmission images of paper (a) 30 min.after deposition of a 50 wt% TEG droplet ( = 5 μl), and (b) after deposition and subsequent evaporation of two water droplets ( = 1 μl) in its center [21].(c,d) Corresponding numerical simulations with grayscale representing the co-solvent content in the pores.Image widths 12 mm.The dashed line in (d) indicates the initial distribution of the added water.

Table 1
Default values of key parameters.initial ink distribution with radius  0 ≤   is present in the pores given by