Superfast non-linear diffusion: Capillary transport in particulate porous media.
The migration of liquids in porous media, such as sand, has been commonly considered at high saturation levels with liquid pathways at pore dimensions. In this letter we reveal a low saturation regime observed in our experiments with droplets of extremely low volatility liquids deposited on sand. In this regime the liquid is mostly found within the grain surface roughness and in the capillary bridges formed at the contacts between the grains. The bridges act as variable-volume reservoirs and the flow is driven by the capillary pressure arising at the wetting front according to the roughness length scales. We propose that this migration (spreading) is the result of interplay between the bridge volume adjustment to this pressure distribution and viscous losses of a creeping flow within the roughness. The net macroscopic result is a special case of non-linear diffusion described by a superfast diffusion equation (SFDE) for saturation with distinctive mathematical character. We obtain solutions to a moving boundary problem defined by SFDE that robustly convey a time power law of spreading as seen in our experiments.
Extremely low volatility (persistent) liquids can spread significantly in porous substrates and wet large volumes over long periods of time before they are removed by evaporation. For example, we found in experiments that a drop of Tris(2-ethylhexyl) phosphate (TEHP - vapour pressure ) deposited on a bed of Ottawa sand, covered a volume of in days, where the spreading essentially came to a stop. With a measured (wet-bed) porosity of , this suggests an average saturation level of ( is the liquid volume divided by the volume of available pore space) and, if uniformly distributed on the surface of spherical grains, a ”film” thickness of . Clearly we are confronted with the mechanics of a process unlike anything found in previous investigations of liquid spreading in environmental porous media.
Equilibrium considerations at low saturation levels have been made in connection with structural properties of granular materials (sands, soils), and these provide a useful starting point [1–3]. The experiments conducted with spherical glass particles (, unknown roughness), shaken with water and allowed to settle, have revealed that liquid domains may exhibit different morphologies. At the low end of saturations the liquid forms isolated pendular rings (liquid bridges) at the points of particle contact. The minimal saturation at the onset of bridge formation was found to be , which is also the level at which the system loses its cohesive property. At higher saturations, , some bridges coallesce into more complex structures (trimers, pentamers and heptamers), but still global connectivity at the pore scale is lacking. [We have observed similar behavior with sand particles at (Fig. 1).]
We anticipate that these overall guidelines for loss of pore-scale connectivity also apply to non-equilibrium systems, such as in our experiment described above, thereby deducing that any transport observed at such low levels of saturation must be due to capillary action at the surface-roughness scale of individual grains (a web-like surface flow geometry). This idea is pursued here by means of a relatively simple mathematical model that ties together the key physics of this transport phenomenon and by highly sensitive experiments devised specifically for accurate tracking of wetting fronts at such trace-quantity conditions.
Previous work with liquids spreading along microscopic surface grooves of various shapes and dimensions has shown that the flow obeys a local Darcy-like law with the permeability related to the groove geometry . Taking this as an approximate model of the micro-flow within the surface roughness, applying intrinsic liquid averaging and the spatial averaging theorem , one can deduce that the macroscopic, averaged over the effective continuum liquid flux and pressure also obey a Darcy-like relationship
Here, and are the liquid velocity and pressure averaged within a surface “groove”, is viscosity, is the area of entrances and exits for the liquid phase through the macroscopic surface element , and is the individual-groove permeability. The latter, if modified by the surface ratio , yields an effective porous medium permeability . An order of magnitude estimate of this quantity can be obtained from the geometry, roughness amplitude and the assumption that the wetted grooves are completely filled [4,6], namely, , where is the effective particle radius. Similarly, the saturation due to the liquid within the roughness is .
The connection of (1) to the total saturation is made by accounting for the liquid content in the bridges . Based on the results of , we obtain, for the case of two identical spheres in contact and assuming near-complete wetting systems, an approximate expression with maximum error for the normalized mean curvature as a function of the normalized bridge volume , with , and . The saturation due to the liquid in the bridges is , and is the particle coordination number. Since , one obtains
where is surface tension.
where and are the characteristic saturation and the length scale of the system and is the diffusion coefficient.
It can be readily seen that the character of this diffusion problem depends utterly on the type of non-linearity dictated by the value of . If , we have what is known as the Porous Medium Equation (PME)  - the typical modelling paradigm in the area of porous transport. Diffusivity in the case of PME yields zero flux at the edge of any initial condition with compact support. This in turn leads to the so-called waiting times and infinite slopes in solutions at moving boundaries, which require special numerical treatment based on a moving-mesh method . If, on the other hand, , which is the case here, the diffusivity, and therefore the flux, become arbitrarily large as ; this is a manifestation of the capillary pressure singularity in (2). In this case, one must have non-zero saturation at the front for a solution to a boundary problem to exist. Accordingly, we must supplement (3) with the front velocity, which from the flux continuity at the moving front with a normal vector , in the dimensionless form reads
A mathematical discussion of this special class of problems can be found in  under the term superfast diffusion (SFD). In our case parameter is fixed and any solution to (3)-(4) belongs to a two-parameter family of functions depending on and . Since is the terminal saturation, its value expresses the final extent of the wet volume. The second parameter is the measure of capillary suction at the front - the larger , the smaller the capillary pressure. The structure of this solution space is illustrated after the presentation of the experimental data, further below.
Our experiments covered a variety of sands and several well-wetting, low-volatility (organophosphate) liquids of varying viscosities and surface tensions ( Tributyl phosphate (TBP), CAS 126-73-8; TEHP, CAS 78-42-2; and Tricresyl Phosphate (TCP), CAS 1330-78-5; Sigma-Aldrich). The liquids have viscosities at , surface tensions, measured in our laboratory at , and contact angles measured on smooth/rough glass respectively. Initial tests showed that spreading is slowed down in the presence of re-entrant microscopic cavities on the sand grains as found for example in erosion sands - they act as irreversible liquid reservoirs. Therefore, it is sufficient for our purposes here to focus on the well-characterized (for example, ) Standard Ottawa sand (EMD Chemicals, p/n SX0075). The sand has an average grain size of and surface roughness in the range . Natural packing (slight shaking to level the bed) yields porosities of , but when wetted we measured (by analysis of images such as in Fig.1) , and a coordination number . The experiments were carried out by depositing gently (fall distance ) a liquid drop and following the evolution of the wet region by imaging the top surface of the bed. In addition, penetration into the bed was observed in half-symmery tests as illustrated in Fig. 2. The time-lapse photography was computer-controlled in a setting that was left undisturbed inside a dark-curtain envelope over periods of many days at controlled temperature. Variability was checked by testing two side-by-side drops placed far enough to avoid interference during spreading, as well as by repeated tests with single drops of two different sizes.
The crucial part of testing was ensuring that the position of the spreading front was captured reliably both in the images and in their analysis. Proper visualization was achieved by UV-excited fluorescence; Coumarin dye, liquid properties confirmed to be unaffected. Photographs were taken by two MP b/w digital cameras (Lumenera Corporation) equipped with macrolenses and focused to resolve individual grains.The lenses were fitted with longpass glass filters to cut off scattered excitation light. No significant background signal could be detected in the absence of the dye in the range of exposures used in the experiments. We have taken images at different values of the exposure to ensure confident capture of the wetted front position at saturation levels down to (mixed-to-equilibrium sand). In analysis, a threshold level of signal at the wetted front was clear and used in automatic scanning of images to determine wetted areas as a function of time.
The results of seven tests are summarized in Fig. 3. All were taken at full symmetry, with various liquids and amounts as quoted in the figure. Inside the SFD regime (marked on the figure) all data fall neatly on a power law with exponent . The small amount of scatter within each run suggests that the slight parallel shifts seen between the individual runs are due to small systematic errors such as deviations in the void fraction and liquid mass values used in the normalization. An estimated error in both these values would result in the error bar shown in Fig. 3. As expected, the SFD regime is attained at , while deviations from the power law set in at a saturation of , and a limiting value of is reached asymptotically over a period of many days. This may be due to the gradual decrease of active contacts between the grains (failure to form bridges due to insufficient liquid), which effectively diminishes the area ratio in (1). For an entire run signal intensity is rather uniform with individual grains appearing illuminated as hypothesized in the theoretical development. Additional tests carried out at half-symmetry have confirmed that in the SFD regime wetted volumes were hemispherical (Fig. 3, inset) and followed the same time law.
To compare experimental results with what is predicted by our SFD model, problem (3)-(4) has been solved numerically at assuming spherical symmetry by a semi-implicit moving-mesh method . The basic character of the solutions is illustrated in Fig 4 (inset). We start simulations at the entrance to the SFD regime with a nearly flat initial saturation profile (, ); the initial wet area provides the length scale and . At high/low capillary suction (low/high ) the spreading goes faster/slower at early times, but in all cases it develops into the -power law seen in the experiments (). The effect of large/small value of is to move down/up the final wet volume. The experimental data for the TEHP drops were placed in Fig. 4 by measuring time from the moment when and scaling it with calculated at and . The numerical solution has been obtained by selecting , according to the experimental data, and . The unique combination of and has been found by variations of parameter . This is possible since both and are monotonic increasing functions of and, as one can see from Fig. 4 (inset), they both result in a parallel shift of the evolution curve in the same direction (the smaller/larger , the smaller/larger both and , so the shift is to the right/left). The effective roughness amplitude (”film” thickness) for the TEHP experiments was found to be (with ). Another case chosen for comparison is the evolution of an TCP drop, which as may be noted stands alone in Fig. 3. The data and numerical solution are shown in the inset to Fig. 4. The TCP case requires and somewhat larger (smaller capillary pressure at the front) corresponding to (with ). This is consistent with the higher values of the contact angle found for this liquid. The obtained non-dimensional characteristic time scale , Fig. 4, suggests that the time should be actually scaled by rather than by . That is substantiating the role of parameter and the capillary pressure in the process of spreading. The obtained values of are consistent with the effective ”film” thickness estimate and suggest that the capillary action occurs in the finest range of the scales found in the roughness of the sand grains .
In conclusion, we have found that highly persistent liquids deposited on porous/granular media spread out to extremely low saturations in a transport regime not previously observable with normal liquids. The process is governed by an equation that belongs to the superfast diffusion class, and the result is a power law for the wetted volume as found in the experiments. The key physics is capillary-driven creeping flow within the roughness-defined channels, which are connected at the grain contacts by liquid bridges. The spreading stops at saturation , which is close to the level at which wet sand loses structural integrity. This suggests the crucial role of the bridges. The key parameter of the problem, , expresses capillary suction at the wetting front and it is found to be related to the contact angle (surface energies).
Acknowledgement. This work is supported by the Joint Science and Technology Office, Defense Threat Reduction Agency (JSTO/DTRA), Threat Agent Science (TAS) of the US Department of Defense. Special thanks to Dr. C.-H. Chang and Dr. S. Cai for the processing of MicroXCT data.
 M Scheel et al. (2008) Nature Materials 7 189–193.
 M Scheel et al. (2008) J. Phys. Condens. Matter 20 494236.
 S. Herminghaus (2005) Advances in Physics 54 221-261.
 Rye, R.R., Yost, F.G., OToole, E.J., (1998) Langmuir 14 3937–3943.
 Whitaker, S., (1969) Ind. Eng. Chem. 61 14-28.
 Bico, J., Thiele, U., Quere, D., (2002) Colloid Surface A 206 41-46.
 Orr, F.M., Scriven, L.E., Rivas, A.P., (1975) J. Fluid Mech. 67 723–742.
 J.L. Vazquez, Smoothing and Decay Estimates for Nonlinear Diffusion Equations - Equations of Porous Medium Type Oxford University Press (2006)
 M.J.Baines, M.E.Hubbard, P.K.Jimack, (2011) Comm. Comput. Phys. 10 509–576..
 The Center for Research Information, Inc., 9300 Brookville Rd, Silver Spring, MD 20910, Health Effects of Trioctyl Phosphate Contract No. IOM-2794-04-001 of the National Academies (2004). D. G. Tuck, Trans. Faraday Soc. 57, 1297, (1961). R. D. Deegan, R. L. Leheny, N. Menon, S. R. Nagel, and D. C. Venerus, J. Phys. Chem. B 103, 4066, (1999).
 Alshibli, K.A., Alsaleh, M.I., (2004) J. Comput. Civil Eng. 18 36–45.