Micromechanics modeling the solute diffusivity of unsaturated granular materials

International Journal of Multiphase Flow 79 (2016) 1–9

Contents lists available at ScienceDirect

International Journal of Multiphase Flow journal homepage: www.elsevier.com/locate/ijmultiphaseflow

Micromechanics modeling the solute diffusivity of unsaturated granular materials Rongwei Yang a,c,∗, Kefei Li a, Eric Lemarchand b, Teddy Fen-Chong c,d a

Civil Engineering Department, Tsinghua University, 100084 Beijing, PR China Université Paris-Est, Laboratoire Navier (UMR 8205), CNRS, Ecole des Ponts ParisTech, IFSTTAR, Champs-sur-Marne 77455, France c Université Paris-Est, Laboratoire Navier (UMR 8205), CNRS, Ecole des Ponts ParisTech, IFSTTAR, Marne-la-Vallée 77420, France d Université Paris-Est, MAST, FM2D, IFSTTAR, Marne-la-Vallée 77447, France b

a r t i c l e

i n f o

Article history: Received 13 May 2015 Revised 1 October 2015 Accepted 21 October 2015 Available online 28 October 2015 Keywords: Unsaturated Solute diffusion Micromechanics Granular material Water film Connectivity

a b s t r a c t This work is devoted to modeling the evolution of the homogenized solute diffusion coefficient within unsaturated granular materials by means of micromechanics approach. On the basis of its distinct role in solute diffusion, the liquid water within unsaturated granular materials is distinguished into four types, namely intergranular layer (interconnected capillary water), isolated capillary water, wetting layer and water film. Application on two sands shows the capability of the model to accurately reproduce the experimental results. When saturation degree is higher than the residual saturation degree Srr , the evolution of homogenized solute diffusion coefficient with respect to the saturation degree depends significantly on the connectivity of the capillary water. Below Srr , depending on the connectivity of the wetting layer, the homogenized solute diffusion coefficient within unsaturated sands decreases by 2–6 orders of magnitude with respect to that in bulk liquid water. The upper bound of the solute diffusion coefficient contributed by the water films is 4–6 orders of magnitude lower than that in bulk liquid water. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction Solute diffusion is an essential transport mechanism in either saturated porous media with low permeability or unsaturated porous media at low saturation degree (Hu and Wang, 2003). Solute diffusion in unsaturated granular materials plays a significant role in agriculture engineering, contaminant transport and remediation, and nuclear waste disposal (Romkens and Bruce, 1964; Conca and Wright, 1992; Hu and Wang, 2003). The liquid distribution within unsaturated granular material is believed to greatly influence the transport properties of the unsaturated granular materials. Solute diffusion at different saturation regimes, namely high saturation regime, intermediate saturation regime and low saturation regime, exhibits distinct transport mechanisms (Hu and Wang, 2003). Therefore, characterizing the evolution of liquid distribution at various saturation degrees is of prior importance in transport problems. Homogenized transport properties are function of liquid content (or saturation degree), which is related to the matric suction (or capillary pressure) according to the water retention curve. In the past decades, extensive empirical models were developed to estimate

∗ Corresponding author at: Civil Engineering Department, Tsinghua University, 100084 Beijing, PR China. Tel./fax: +86 10 6279 7993. E-mail address: [email protected] (R. Yang).

http://dx.doi.org/10.1016/j.ijmultiphaseflow.2015.10.004 0301-9322/© 2015 Elsevier Ltd. All rights reserved.

the effective solute diffusion coefficient within unsaturated porous media (Romkens and Bruce, 1964; Mehta et al., 1995; Olesen et al., 1999; Hu and Wang, 2003). These models are phenomenological since local physical phenomena are rarely taken into account. For instance, below residual saturation degree, the thick water films (the wetting layer) attached on the solid grain surface are found to be dominant in the transport problems, which is the main physical phenomenon neglected by the Archie’s law (Conca and Wright, 1992; Hu and Wang, 2003; Han et al., 2009). Therefore, micromechanics approach, accounting for both microstructure morphology and local phase properties, has been used widely in estimating the homogenized transport properties of porous media (Gruescu et al., 2007; Lemarchand et al., 2009; Nguyen, 2014). The latter often resorts to the classic Eshelby-based schemes such as self-consistent scheme (SC) and Mori-Tanaka scheme (MT) (Lemarchand et al., 2009). When MT is used in the transport problems, the capillary water is usually considered as matrix, the solid phase inclusions along with gaseous phase inclusions are embedded in the matrix. This morphology ensures the interconnectivity of the capillary water during the overall desaturation (drainage) process, which is in contradiction with several experimental results which exhibit notable percolation threshold at low saturation degree (Romkens and Bruce, 1964; Hamamoto et al., 2010). As an alternative, SC is used to estimate the homogenized solute diffusion since the percolation threshold can

2

R. Yang et al. / International Journal of Multiphase Flow 79 (2016) 1–9

Fig. 1. Schematic illustration of liquid distribution within unsaturated granular material at distinct saturation regimes.

be inherently accounted for in this scheme. Unfortunately, within unsaturated porous media, the percolation threshold saturation degree Srp estimated by SC is 31φ (φ is pore volume fraction and φ ≥ 1/3)

(Yang, 2013), which is obviously overestimated (Koelman and De Kuijper, 1997). The overestimation lies in the fact that the connectivity of the capillary water within unsaturated porous media, as an important topological information, is not physically incorporated into the Eshelby-based micromechanics estimates (Kleinfelter-Domelle and Cushman, 2012). In this study, four different types of liquid water, namely intergranular layer (interconnected capillary water), isolated capillary water, wetting layer and water film, are distinguished and specified within unsaturated granular materials. Based on the local liquid characterization, a physical-based micromechanics model is developed to predict the solute diffusion coefficient of the unsaturated porous media. The micromechanics model is validated by the experimental results from the literature and compared with the empirical Archie’s law. The model provides a link between the local parameters such as connectivity of capillary water and different tendencies of the homogenized solute diffusion coefficient of the two sands with respect to saturation degree.

2. Physical characterization of unsaturated granular materials 2.1. Saturation regimes On the premise of the two critical saturation degrees Srr and 1%, three saturation regimes in unsaturated porous media, namely high saturation regime, intermediate saturation regime and low saturation regime, are used to distinguish the liquid configuration at different saturation degrees. At high saturation regime (Sr ≥ Srr ), the capillary water can be classified as percolating (interconnected) capillary water and isolated capillary water (Raoof and Hassanizadeh, 2013). The latter is attributed to either the segmentation by the percolating non-wetting gaseous phase (Raoof and Hassanizadeh, 2013) or the existence of the dead-end pore space (Nakashima, 1995). Therefore, at high saturation regime, the capillary water is classified as isolated capillary water and interconnected capillary water. The latter is believed to play a vital role in the solute diffusion and is idealized as an intergranular layer attached on the surface of solid grain in our model (see “High saturation regime” in Fig. 1).

Intermediate saturation regime corresponds to a regime where the saturation degree is below the residual saturation degree Srr (Sr < Srr ) (Romkens and Bruce, 1964; Fredlund et al., 1997). At intermediate saturation degree, the intergranular layers disappear since the connectivity of the capillary water vanishes. As shown in Fig. 1, the capillary water trapped in pendular rings and surface roughness forms a continuous transport pathway, which will play an important role in the solute diffusion (Or and Tuller, 2000; Blunt et al., 2002). This continuous transport pathway (the thick water film) is defined as “wetting layer” by Blunt et al. (2002), which keeps interconnected even when the saturation degree Sr reaches as low as 1% (Dullien et al., 1989). The thickness of the wetting layer is closely associated with the morphological characteristics of surface roughness and the capillary pressure. Owing to the complicated morphological characteristics of the surface roughness, it is a difficult task to precisely characterize the relationship between the capillary pressure and thickness of the wetting layer (Or and Tuller, 2000; Blunt et al., 2002; Kibbey, 2013). However, it is generally agreed that the magnitude of the thickness of the wetting layer is the order of the characteristics size of surface roughness, ranging between microns to sub-microns (Dullien et al., 1989; Han et al., 2009; Kibbey, 2013). In addition, there still exists certain amount of isolated capillary water (Kibbey, 2013), either in small capillary pores or in pendular rings (or surface pits and concaves), which is bridged by the wetting layers (see “intermediate saturation regime” in Fig. 1). As illustrated in “Low saturation regime” in Fig. 1, with the progressive drying (e.x. Sr < 1% (Dullien et al., 1989; Blunt et al., 2002)), the low saturation degree regime is reached: capillary water trapped in small capillary pores (the isolated capillary pore water) is totally drained out. Correspondingly, the interconnected wetting layers disappear. Thereafter, the capillary water trapped in surface roughness and pendular rings is connected by the water films. The latter predominate the solute diffusion at low saturation regime. In the following context, the superscripts/subscripts s, g, il, ip, wl and f denote solid phase, gaseous phase, intergranular layer, isolated capillary water, wetting layer and water film, respectively. The boldfaced letters (e.g. A) stand for second order tensor while the underlined letters (e.g. H ) denote vectors. 2.2. Thickness of water film At nano-scale, the stability and thickness of water film are determined by the so-called disjoining pressure (h) (Derjaguin and Churaev, 1978; Israelachvili, 1991). The latter is linear sum of three different components: the relative long range repulsive electrostatic force e (h), the Van der Waals component v (h), the structural component of disjoining pressure s (h) (Churaev and Derjaguin, 1985; Majumdar and Mezic, 1999; Gonçalvès et al., 2010):

(h) = e (h) + v (h) + s (h)

(1)

with (Churaev and Derjaguin, 1985; Israelachvili, 1991):

⎧ ⎨e (h) = 64nkT ζ1 ζ2 exp( − κ h) A  (h) = − 6π h3 ⎩ v s (h) = Ksr exp( − h/λsr ) + Klr exp( − h/λlr )

(2)

where Debye length κ1 = (0 kT )1/2 (2z2 e2 n)−1/2 ;  (=80) is the dielectric permittivity of liquid water,  0 (=8.854 ×10−12 C2 J−1 m−1 ) is the permittivity of the free space; k (=1.381 ×10−23 J · K−1 ) is the Boltzmann constant; T is the temperature (in Kelvin); n is the number density of the ions, n = Na × ρ , ρ (in mol.L−1 ) is the concentration of the bulk electrolyte solution, Na (=6.022 × 1023 ) is Avogadro’s constant; e (=1.602 × 10−19 C) is the elementary charge; z is the electron charge; ζi =

exp(zeψi /(2kT ))−1 exp(zeψi /(2kT ))+1

(i ∈ {1, 2}), ψ i are the surface poten-

tials of ith interface (i ∈ {1, 2}), herein, subscript 1 denotes solid-water film interface, subscript 2 denotes gas-water film interface; A is the

R. Yang et al. / International Journal of Multiphase Flow 79 (2016) 1–9

3

Fig. 2. Morphological representation of each phase within unsaturated granular materials when Sr > Srr (a) or Sr < Srr (b); (b) can be extended to low saturation regime (Sr < 1%) when the interconnected wetting layers (in carmine) vanish; the problems can be separated into three Eshelby-based problems as shown in the insets; the yellow color denotes the homogenized medium we are looking for, the other colors represent each phase as shown in Fig. 1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

Hamaker constant; Ksr , Klr , λsr and λlr are the four fitting parameters for structural component of disjoining pressure; h is the thickness of the water film. When uniform thickness of water film on the perfectly wetted solid surface is assumed, disjoining pressure can be either related to matric suction  or related to the relative humidity hr by Kelvin equation (Israelachvili, 1991):

(h) = e (h) + v (h) + s (h) = − = −

2γlg RT ln(hr ) + Vl Rs

(3)

where Vl is the molar volume of the liquid water, R (=8.314 −1 J.mol .K−1 ) is the perfect gas constant. Substituting Eq. (2) into Eq. (3) yields an implicit equation, according to which the thickness of water film adsorbed on the solid grains can be estimated. 2.3. Local solute diffusion within water film Local solute diffusion through porous media exhibits multiphysics behaviors: the solute diffusion coefficient within capillary water can be treated as that in bulk liquid water; while the reduced solute diffusivity within thin water film was observed (Grathwohl, 1998; Bourg et al., 2007). The reasons for the reduced solute diffusivity may lie in the increased viscosity of water film, the ionic adsorption and interaction with the charged solid surface (Conca and Wright, 1992). Generally, a scaling constrictive factor δ for solute τ is introduced to characterize the reduced solute diffusivity within water film with respect to Dτ (Bourg et al., 2007):

Dτf =

δ Dτ

with

δ≤1

(4)

where Dτf , Dτ are the diffusion coefficients of solute species τ in water film and bulk liquid water, respectively. In order to determine the constrictive factor δ , detailed knowledge about the dependence of the viscosity on thickness of water film, ionic interaction, mineralogy and surface properties of the granular material is required. However, it can be assumed that δ = 1 when h ≥ 10 nm, while when h < 10 nm, δ < 1 (Grathwohl, 1998). In this study, δ = 1 is adopted in the following model for simplification. 3. Theoretical framework According to Section 2.1, part of the “liquid water” is always attached to the pore/solid interface during the desaturation process.

This “liquid water” varies from intergranular layer at high saturation degree, wetting layer at intermediate saturation degree to water film at low saturation degree. Hence, as depicted in the insets of Fig. 2, intergranular layer + wetting layer + solid grain composite and wetting layer + water film + solid grain morphologies are sufficient to represent the morphologies of liquid water attached on the solid grain during the desaturation process. Based on the characterization of liquid distribution within unsaturated granular material in Section 2.1, a schematic illustration of the representative elementary volume (REV) of unsaturated granular material (with domain ) is shown in Fig. 2. The REV is composed of the solid phase (domain s ), the gaseous phase (domain g ), the isolated capillary water (domain ip ), the wetting layer (domain wl ), the intergranular layer (domain il ) or the water film (domain f ). That is to say, when Sr > Srr ,  = s ∪ g ∪ ip ∪ il ∪ wl ; while when Sr < Srr ,  = s ∪ g ∪ wl ∪  f . A uniform macroscopic concentration gradient H is prescribed on the boundary of the REV (∂ ) in the sense of Hashin (Dormieux et al., 2006). The average of a local physical quantity a( z ) over  can be expressed as (Dormieux et al., 2006):

a=



1

||



a(z)d

(5)

The local diffusion problem in this unsaturated granular material can also be defined over a REV (domain ) by the following set of equations (Dormieux et al., 2006):



τ

div j = 0 τ j = −D(z) gradρ τ ρτ = H · z

with:

D(z) =



Ds = Dg = 0 Dτ Dτf = δ Dτ

z∈ z∈ z ∈ ∂

(6)

∀z ∈ s ∪ g ∀z ∈ ip ∪ wl ∪ il ∀z ∈  f

(7)

where j τ is the solute mass flux, z is the position vector, D( z ) is the local solute diffusion coefficient, ρ τ is the concentration of the solute species τ and grad ρ τ is the concentration gradient of the solute species τ . The concept of second order concentration tensor A( z ) is introduced in order to express the heterogeneous solution of Eq. (6) as

4

R. Yang et al. / International Journal of Multiphase Flow 79 (2016) 1–9 Table 1 Characteristics of the two sands (Romkens and Bruce, 1964; Lim et al., 1998).

(Dormieux et al., 2006):

gradρ τ = A(z) · H

(8)

where the concentration tensor A( z ) is a function of the geometrical microstructure. In terms of the uniform boundary condition (third formula in Eq. (6)) and Eq. (5), we have (Dormieux et al., 2006):

gradz ρ τ = H ⇔ A = 1

(9)

where 1 is the second order unit tensor. According to the upscaling rule defined in Eq. (5), the macroscopic mass flux can be determined as Jτ = jτ (z) (Dormieux et al., 2006). Inserting second equation of Eqs. (6) and (8) into Eq. (5) yields a macroscopic Fick’s law (Dormieux et al., 2006):

Jτ = jτ (z) = −Dhom .H

with: Dhom = D(z)A

(10)

where Dhom is the homogenized solute diffusion coefficient tensor. If the local and macroscopic isotropies of unsaturated granular material are assumed, we have Dhom = Dhom 1, Ai = 13 1 : Ai , where Ai is the average concentration gradient coefficient of ith phase, i ∈ {s, g, ip, il, wl, f}. As introduced in Sec. 4.1.2 and shown in Fig. 2, when Sr ≥ Srr (at high saturation regime), the liquid within unsaturated granular materials can be separated into the intergranular layer, isolated capillary pore water and wetting layer; similarly, when Sr < Srr (at intermediate saturation regime), the liquid can be distinguished into the wetting layer, isolated capillary pore water and water film. The homogenized solute diffusion coefficient Dhom can be expressed as:

⎧ ⎪ hom ⎪ ⎪ ⎨D = ⎪ ⎪ ⎪ ⎩Dhom =

φil Ail Dτil + φip Aip Dτ + φwl Awl Dτ when Sr ≥ Srr φil Ail + φip Aip + φwl Awl + φs As + φg Ag φ f A f Dτf + φip Aip Dτ + φwl Awl Dτ φ f A f + φip Aip + φwl Awl + φs As + φg Ag

when Sr < Srr (11)

with

 φil Ail + φip Aip + φwl Awl + φs As + φg Ag = 1 when Sr ≥ Srr φ f A f + φip Aip + φwl Awl + φs As + φg Ag = 1 when Sr < Srr

(12)

where φ i is the average concentration coefficient and volume fraction for ith phase (i ∈ {s, g, ip, il, wl, f}). When φwl = 0 at Sr < Srr , the intermediate saturation regime transforms into low saturation regimes (i.e., Sr < 1%) and water films govern the solute diffusion process, which is introduced in Section 2.1. If Ai (i ∈ {s, g, ip, il, wl, f}) is known, Dhom may be exactly determined. However, the exact determination of the microstructure (Ai ) is a complicated task: we rather resort to estimates. The granular materials do not exhibit clear matrix+inclusion morphologies since any phase within these materials not play the role of the embedding matrix. Their apparent polycrystalline nature rather suggests a perfect disordered microstructure, which is classically accounted for by the well-known self-consistent scheme (SC). SC consists of assuming that each particle of a given phase (pore or solid) reacts as if it were embedded in the equivalent homogeneous medium which is looked for (Dormieux et al., 2006). The SC instead of MT is well adapted for the granular materials based on two further reasons: (1) SC is the sole estimate which is able to discuss percolation effects likely to occur within porous media during the desaturation process; (2) the interconnected capillary pore water (intergranular layer) should be taken as matrix in Mori-Tanaka scheme; in this case, the morphology of the isolated capillary pore water can not be accounted for in Mori-Tanaka scheme. In SC, the homogenization of solute diffusion can be decomposed into three auxiliary Eshelby-based problems (see the insets in Fig. 2). When Sr ≥ Srr or Sr < Srr , the average concentration coefficient of

Property

Unit

Beaver Creek sand

Romkens sand

Porosity Dry bulk density Tracer Initial concentration of tracer Residual saturation degree Srr Dτ in bulk liquid water Saturated diffusion coefficient

% Kg · m−3 mol · L−1 % m2 · s−1 m2 · s−1

38.2 1647 ± 7 − Cl 0.02 9 2.03 × 10−9 0.535 × 10−9

42.2 1537 ± 19 NO−1 3 0.1 7 1.90 × 10−9 0.590 × 10−9

each phase can be solved by the Eshelby-based solution: there exists a uniform concentration gradient within the single imbedded spherical inclusion (Eshelby, 1957). The average concentration gradient coefficients of each phase Ai (i ∈ {s, g, ip, il, wl, f}) are given in the Appendix A. For detailed derivation one can refer to Yang (2013). According to the expressions of Ai (i ∈ {s, g, ip, il, wl, f}) listed in the Appendix A, it is readily found that Ai are the functions of volume fractions and local solute diffusivities of each phase. Therefore, the estimate of Dhom turns to the determination of the volume fraction of each phase φ i . 4. Application: nonsorbing solute diffusion in unsaturated sands The micromechanics model for solute diffusion is used to simulate the evolution of the homogenized solute diffusion within unsaturated sands. Two sets of the solute diffusion experiment within unsaturated sand have been carried out by Lim et al. (1998) and Romkens and Bruce (1964). A detailed description of the procedures for diffusion testing is given in their works (Romkens and Bruce, 1964; Lim et al., 1998). To avoid the ions-surface interaction, the ex+ perimental results of nonreactive Cl− and NO−1 3 instead of K are adopted to compare with the model results. Firstly, the main characteristics of two sands are listed in Table 1. In order to manage the granulometry of the sand, we introduce a continuous grain size distribution function in the model. The latter will account, in a first approximation, for the specific microstructure of granular material. Moreover, the water retention curves of the two sands are introduced in the micromechanics model to quantify the volume fractions of capillary water during desaturation process. It is assumed that, during the desaturation process, the capillary pores are emptied from the bigger ones to the smaller ones in a well organized arrangement (Dullien, 1992). 4.1. Phase volume fractions Two empirical models for water retention curves at Sr ≥ Srr (Fredlund and Xing, 1994) and at Sr < Srr (Campbell et al., 1992) are used to characterize the evolution of the saturation degree with the matric suction (or relative humidity):

⎧ ⎛ ⎞m ⎪ ⎪ ⎪ ⎪ ⎜ ⎟ ⎪ 1 ⎪ ⎪ Sr = Srr + (1 − Srr )⎜  n ⎟ ⎪ ⎝ ⎠ ⎪ − ⎪ ⎨ ln[exp(1) + ] ⎪ whenSr ≥ Srr (a) ⎪ ⎪ ⎪ ⎪ log10 ( −  /re f ) ⎪ r ⎪ Sr = Sr 1 − ⎪ ⎪ 5 ⎪ ⎩ when Sr < Srr (b)

k

(13)

where  (in kPa) is the matric suction and is associated with relative humidity in terms of Eq. (3); k, m n and Srr are the fitting parameters for two kinds of sand, which are listed in Table 2;  ref is the reference matric suction taken as re f = 10 kPa (Campbell et al., 1992).

R. Yang et al. / International Journal of Multiphase Flow 79 (2016) 1–9

5

Fig. 3. Water retention curves for Beaver Creek sand, Romkens’ sand and Shonai sand; the matric suctions associated with the residual saturation degrees for Beaver Creek sand (9%), Romkens’ sand (7%) and Shonai sand (13%) are −20 kPa, −30 kPa and −20 kPa, respectively; experimental results are after Romkens and Bruce (1964), Lim et al. (1998) and Mehta et al. (1995).

by Eq. (13)(b), φ f is the volume fraction of water film. In summary, we have:

Table 2 Fitting parameters for two kinds of sand. Parameters

k

n

m

Srr

Beaver Creek sand (Lim et al., 1998) Romkens’ sand (Romkens and Bruce, 1964) Shonai sand (Mehta et al., 1995)

3 5.42 1.934

10 12.06 3.794

2 1.94 3.100

9% 7% 13%

The fitting water retention curves for Beaver Creek sand and Romkens’ sand are depicted in Fig. 3. The model is validated by the experimental results when Sr > Srr . Unfortunately, the experimental data of Beaver Creek sand and Romkens’ sand are unavailable when Sr < Srr . Therefore, the experimental data of Shonai sand (Mehta et al., 1995) are used to validate the Campbell’s model (see the second equation in Eq. (13)) when Sr < Srr . It can be found that the model results agree well with the experimental data of Shonai sand. The volume fraction of the solid grains φ s can be readily determined as φs = 1 − φ , where φ is the volume fraction of the pore space (porosity). Similarly, the volume fraction of the gaseous phase φ g can be easily determined as φ(1 − Sr), Sr is associated with the matric suction (or relative humidity) by Eq. (13). 4.1.1. Intergranular layer, isolated capillary water and wetting layer From the topological point of view, the volume fraction of the intergranular layer (interconnected capillary water) should be determined by experimental devices such as micro-CT or environmental scanning electron microscopy (ESEM) (Wang et al., 2013, 2014). However, no other additional topological information such as the connectivity of the capillary water is given in the Romkens and Bruce (1964) and Lim et al. (1998). An empirical expression β Srχ is thus introduced to quantify the volume fraction of intergranular layer φ il when Sr > Srr (Weiss et al., 2012). β and χ are the empirical parameters which characterize the connectivity of the capillary water during desaturation. Therefore, φil = β Srχ (φ Sr − φwl ) and the volume fraction of the isolated capillary water φip = (1 − β Srχ )(φ Sr − φwl ) are determined, where Sr is related to the capillary pressure by Eq. (13)(a). At the high saturation regime, the thickness of the wetting layer is assumed to be 0.6 μm (Han et al., 2009). Hence, φwl = 0.6 × 10−6 S, where S (in m2 /m3 ) is the specific surface of the granular material. Likewise, when Sr < Srr , φil = 0, a term η is used to quantify the connectivity of wetting layer (Or and Tuller, 2000), namely φwl = η(φ Sr − φ f ) and φip = (1 − η)(φ Sr − φ f ), in which Sr is determined

 φil = β Srχ (φ Sr − φwl ) Sr ≥ Srr : φip = (1 − β Srχ )(φ Sr − φwl ) φwl = 0.6 × 10−6 S  φil = 0 r Sr < Sr : φip = (1 − η)(φ Sr − φ f ) φwl = η(φ Sr − φ f )

(14)

(15)

At high and intermediate regimes, taking advantage of Eq. (3) and Eq. (13), φ g , φ il , φ ip and φ wl can be explicitly determined by the relative humidity. 4.1.2. Water film In order to estimate the volume fraction of water film, the only information we need to know is the specific surface of the granular materials. The latter can be determined by the grain size distribution function ν (D) as:



Sg =

Dmax

Dmin

6 ν(D)dD = 6 D



Dmax

ν(D) D

Dmin

dD

(16)

where Dmin and Dmax are the minimum and maximum diameters of sand. For Beaver Creek sand, Dmin = 10−3 mm and Dmax = 10 mm. The grain size distribution function ν (D) is related to the cumulative curve function Pp (D) by:

ν(D) =

dPp (D) dD

with (Fredlund et al., 1997):

(17)

⎡

⎡



 ⎤7 ⎤

Dr 1 ⎢ ⎢ D ⎥ Pp (D) =    g gm  ⎣1 − ⎣  Dr ⎦ a ln 1 + ln exp(1) + Dm D ln 1 +

⎥ ⎦

(18)

where the fitting parameters for Beaver Creek sand are determined as ga = 0.2485, gm = 4.8109, gn = 1.7015, Dm = 0.0001, Dr = 36.6213 (Fredlund et al., 1997). Inserting Eqs. (17) and (18) into Eq. (16) yields the specific surface for Beaver Creek sand S = 6.78 × 104 m2 /m3 . Unfortunately, the grain size distribution curve for Romkens’ sand is not given in Romkens’ work. However, we know that the grain sizes of Romkens’ sand distributed in a narrow range of 105 μm to 210 μm (Romkens and Bruce, 1964). In this case, the grain size of Romkens’

6

R. Yang et al. / International Journal of Multiphase Flow 79 (2016) 1–9 Table 3 Parameters for the quartz-NaCl aqueous film-air system at T=293 K. Parameters

Value (unit)

Reference

z

1 −50 mV −25 mV −0.87 ×10−20 J 3 × 108 Pa 0.3 nm 2 × 106 Pa 2 nm

(Tokunaga, 2011) (Tokunaga, 2011) (Israelachvili, 1991) (Churaev and Derjaguin, 1985) (Churaev and Derjaguin, 1985) (Churaev and Derjaguin, 1985) (Churaev and Derjaguin, 1985)

ψ1 ψ2 A Ksr

λsr

Klr

λlr

Here we investigate the evolution of the thickness of water film on flat solid surface (Rs → ∞), the capillary effect induced by grains 2γ

Fig. 4. Relationship of the thickness of water film with relative humidity, M = mol.L−1 ; experimental results for water film-fused quartz system are after Sumner et al. (2004), the thickness of mono liquid layer is 2.8 A˚ (Boyarskii et al., 2002).

sand is assumed to be mono disperse and the diameter of Romkens’ sand is assumed to be D = (105 + 210)/2 ≈ 150 μm. Hence, the specific surface of the Romkens’ sand is S ≈ 6(1 − φ)/D = 2.31 × 104 m2 /m3 . The volume fraction of the water film φ f is readily estimated as:

φ f = hS

size ( Rslg ) is thus omitted in Eq. (3). Adopting the parameters listed in Table 3, the evolution of the thickness of water film with the relative humidity hr is depicted in Fig. 4. As can be seen from Fig. 4, the model results (the red and green curves) are in good agreement with the experimental results at relative humidity hr lower than 99%. When hr < 85%, the effect of the concentration of the solution is negligible for thin water film. This means that the electrostatic component of disjoining pressure is negligible for thin water film, the thickness of thin water film is thus mostly determined by the structural component and Van der Waals component of disjoining pressure. However, the influence of the solute concentration is notable when hr > 85%: the higher the concentration (green curve), the higher the thickness of water film. That is because the electrostatic component is a long range force and it plays a significant role in the total disjoining pressure at high relative humidity.

(19)

The thickness of the water film h may be determined by Eq. (3) and thus φ f can be determined by the relative humidity. 4.2. Results and discussion 4.2.1. Thickness of water film Eqs. (1)–(3) are used to estimate the thickness of water film in quartz-water film-air system, here the solution is NaCl electrolyte. The parameters used are given in Table 3.

(a) Beaver Creek sand

4.2.2. Evolution of Dhom /Dτ with Sr (Sr ≥ Srr ) Sand is a typical granular material which has the simplest microstructure in the natural geomaterials. However, as can be found from Fig. 5, the experimental results of two sands exhibit quite different tendencies: Dhom /Dτ in Beaver Creek sand increases gently with Sr, there is no obvious percolation threshold; on the contrary, the experimental results of Romkens sand exhibit strong percolation threshold effect: Dhom /Dτ in Romkens sand increases rapidly at high saturation degree (Sr ≥ 70%), when Sr < 70%, Dhom /Dτ increases slightly with Sr.

(b) Romkens sand

Fig. 5. Evolution of Dhom /Dτ with Sr when Sr ≥ Srr , twl is the thickness of the wetting layer, experimental results are after Romkens and Bruce (1964) and Lim et al. (1998).

R. Yang et al. / International Journal of Multiphase Flow 79 (2016) 1–9

(a) Beaver Creek sand

7

(b) Romkens sand

Fig. 6. Evolution of Dhom /Dτ with Sr at various η when Sr < Srr , the cyan dash dotted curves are the solute diffusion contributed by the water films. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article).

In our model, the term β characterizes the connectivity of the intergranular liquid at fully saturated case (Sr = 1). As the saturated homogenized solute diffusion coefficients for the two sands are given in Table 1, making Sr = 1 and inserting Eq. (14) into Eq. (11) yields β = 0.606 for Beaver Creek sand (φ = 0.382) and β = 0.727 for Romkens sand (φ = 0.422). As expected, β increases with φ , which means that the connectivity of the intergranular liquid in fully saturated sands increases with φ . The exponent parameter χ is used to quantify the evolution of φ il during desaturation process (see Eq. (14)). The greater the χ , the greater loss of the connectivity of the capillary water (intergranular layer) and thus the more pronounced percolation effect during desaturation. As depicted in Fig. 5, when χ = 1 (for Beaver Creek sand) and χ = 10 (for Romkens sand) as well as twl = 0.6 μm for both sands, the model results (solid curves) agree well with the experimental results of the two sands, respectively. According to the values of χ , it can be inferred the intergranular liquid within Beaver Creek sand is well connected when Sr ≥ Srr while the intergranular liquid within Romkens sand loses its connectivity sharply at high saturation degree (Sr > 70%). As illustrated in Fig. 5, the quite different experimental results of the two sands during desaturation are attributed to the distinct connectivity of intergranular layer (interconnected capillary water) during desaturation process. Further light should be shed on understanding the physical meaning of χ and β . Sophisticated experimental devices such as 3D micro-CT could be applied to better quantify β and η within unsaturated granular materials; in parallel, theoretical tools such as percolation theory can also be used to interpret the physical meaning of these two parameters. Archie’s law, in the form of Dhom /Dτ = φ α Srν (Archie, 1942), is used to compare with the present model results. The fitting parameters of Archie’s law for Beaver Creek sand and Romkens sand are α = 1.38 and ν = 2 and α = 1.36 and ν = 6, respectively. With the choice of χ = 1, twl = 0.6 μm for Beaver Creek sand, and of χ = 10 and twl = 0.6 μm for Romkens sand, the present model results are comparable with the empirical Archie’s law. However, the proposed micromechanics model allows taking account of the several mechanisms introduced in the previous section, i.e., the separation of the liquid water morphologies and the local transport properties, which

can be neglected or not depending on the problem under consideration. As a comparison, the Archies law is essentially empirical and it is not valid below intermediate saturation degree Srr (Han et al., 2009). 4.2.3. Evolution of Dhom /Dτ with Sr (Sr < Srr ) Owing to the detection limit of measuring device as well as being time-consuming and laborious (Conca and Wright, 1992), there are rather scarce experimental results for solute diffusion coefficient below intermediate saturation degree (Sr < Srr ). Herein, we attempt to predict the solute diffusion coefficient below intermediate saturation degree (Sr < Srr ) by means of the micromechanics model. When Sr < Srr and the matric suction reaches as low as −50 kPa, according to Eq. (3), the thickness of the water film is lower than 10 nm. Therefore, the constrictive factor δ < 1 should be accounted for in the model. Nevertheless, for the simplificity, δ is also assumed to be 1 in the model. In this case, it can be inferred that the solute diffusion contributed by the water films is the upper bound. When Sr < Srr , the parameter η is introduced to quantify the connectivity of the wetting layers (see Section 4.1.1). As plotted in Fig. 6, the influence of saturation degree on Dhom /Dτ for the two sands significantly depends on η. The evolution of Dhom /Dτ with respect to Sr increases with η. Depending on the connectivity of the wetting layers (0 < η < 1), Dhom /Dτ (the solid curves in Fig. 6) decreases by 2–6 order of magnitude when Sr < Srr , which is in line with the results of Nakashima (1995) and Koelemeijer et al. (2012). Moreover, the solute diffusion contributed by the water films for the two sands is depicted as cyan dashdotted curves in Fig. 6. It can be found that the upper bond of the solute diffusion coefficient contributed by the water film is 4–6 order lower than that of the bulk liquid water. The impact of the water film becomes significant when the wetting layers are poorly interconnected (e.g. η < 0.05). It should be noted that the reason for the negligible volume effect of the water film is attributed to the low specific surface area of the sand. For other polycrystalline microstructures with high specific surface area (e.g. clay particles arrangement, C-S-H particle arrangement), the water films are anticipated to play an essential role in the solute diffusion. At low saturation regime (Sr < 1%), the interconnected wetting layers vanish and the water films govern the solute diffusion (see Fig. 6).

8

R. Yang et al. / International Journal of Multiphase Flow 79 (2016) 1–9

5. Concluding remarks In this study, the liquid water within unsaturated granular materials is distinguished into four types: intergranular layer (interconnected capillary water), isolated capillary water, wetting layer and water film. On the basis of local characterization of these phases, a micromechanics model is developed to estimate the homogenized solute diffusion coefficient within unsaturated granular materials. Two groups of simulation on the solute diffusion in unsaturated sands are carried out to validate the micromechanics model. When saturation degree Sr is higher than residual saturation degree Srr , two empirical parameters β and χ are introduced to characterize the volume fraction of the intergranular layer (connectivity of the capillary water) during desaturation. β increases with pore volume fraction φ while the percolation effect becomes more and more notable with increasing χ . With given β and χ for the two sands, the model results fit well with the experimental results and are comparable with Archie’s law. The proposed micromechanics model is physically based and it allows taking account of the several mechanisms such as liquid separation and local solute diffusion. This makes the micromechanics model be a alternative tool to estimate the solute diffusion coefficient at intermediate saturation regime and low saturation regime (Sr < Srr ). A parameter η is introduced to quantify the connectivity of the wetting layers. Depending on the connectivity of the wetting layers (0 < η < 1), the normalized solute diffusion coefficient Dhom /Dτ decreases by 2–6 order of magnitude when Sr < Srr . Moreover, the upper bound of the solute diffusion coefficient contributed by the water film is 4–6 order lower than that of the bulk liquid water. However, the former will become significant when the wetting layer is poorly interconnected. This is a promising result since experimental tests on solute diffusion coefficient at intermediate saturation degree (Sr < Srr ) are time consuming.

Appendix A When Sr ≥ Srr , there exist the solid phase (denoted by subscript s), intergranular layer (denoted by subscript il), isolated capillary pore water (denoted by subscript ip), wetting layer (denoted by subscript wl) and gaseous phase (denoted by subscript g). Therefore, when Sr > Srr , the average concentration gradient coefficients of solid phase, wetting layer, and intergranular layer are given in Eq. (A.1):

⎧ (φs + φwl )(φs + φwl + φil )Dhom 27 ⎪ As = ⎪ ⎪ hom 2 15D φs φwl + 3Dτ δφs φwl + Dτ δφil φwl + ... ⎪ ⎪ ⎪ ⎪ ⎪ (φs + φwl )(φs + φwl + φil )Dhom ⎪ ⎪ ⎪ 2 + 3Dτ φ φ + 2Dτ φ φ hom ⎪ ... + 9D φs2 + 6Dhom φwl il s il wl + ... ⎪ ⎨ hom r (φ + φ )(φ + φ + φ ) D s s wl wl il Sr ≥ Sr : ⎪ ⎪ ... + 2Dhom δφil φwl + 6Dhom φil φs + 4Dhom φil φwl ⎪ ⎪ ⎪ ⎪A = 2 As ⎪ wl ⎪ ⎪ 3 ⎪ ⎪ ⎪ 2 δφwl + 3φs + 2φwl ⎪ ⎩Ail = As 9 φs + φwl (A.1) Likewise, when Sr < Srr , the interconnected capillary pore water drained out, there exist the wetting layer, water films (denoted by subscript f), isolated capillary pore water, solid phase and gaseous phase. The average concentration gradient coefficients of solid phase, water film, and wetting layer can be formulated in Eq. (A.2) as:

⎧ (φs + φwl )(φs + φ f + φil )Dhom 27 ⎪ ⎪ As = ⎪ hom ⎪ 2 15D φs φ f + 3Dτ δφs φ f + Dτ δφwl φ f + ... ⎪ ⎪ ⎪ ⎪ (φs + φ f )(φs + φ f + φwl )Dhom ⎪ ⎪ ⎪ ⎪ hom ⎪ ... + 9D φs2 + 6Dhom φ 2f + 3Dτ φwl φs + 2Dτ φwl φ f + ... ⎪ ⎨ (φs + φ f )(φs + φ f + φwl )Dhom Sr ≥ Srr : ⎪ hom ⎪ ... + 2D δφwl φ f + 6Dhom φwl φs + 4Dhom φwl φ f ⎪ ⎪ ⎪ ⎪ 2 ⎪ ⎪ A f = As ⎪ ⎪ 3 ⎪ ⎪ ⎪ δφ + 3φs + 2φ f ⎪ ⎩Awl = 2 f As 9 φs + φ f (A.2) The average concentration gradient coefficients of gaseous phase (resp. isolated capillary pore water) at Sr > Srr and Sr < Srr are the same, which are given in Eq. (A.3):

Aip =

3Dhom Dτ + 2Dhom

Ag =

3 2

(A.3)

References Archie, G., 1942. The electrical resistivity log as an aid in determining some reservoir characteristics. Trans. Am. Inst. Min. Metall. Pet. Eng. 146, 54–62. Blunt, M.J., Jackson, M.D., Piri, M., Valvatne, P.H., 2002. Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow. Adv. Water Resour. 25, 1069–1089. Bourg, I., Sposito, G., Bourg, A., 2007. Modeling cation diffusion in compacted watersaturated sodium bentonite at low ionic strength. Environ. Sci. Technol. 41, 8118– 8122. Boyarskii, D., Tikhonov, V., Komarova, N.Y., 2002. Model of dielectric constant of bound water in soil for applications of microwave remote sensing. Prog. Electromagn. Res. 35, 251–269. Campbell, G., Shiozawa, S., van Genuchten, M.T., Leij, F.J., Lund, L.J., 1992. Prediction of hydraulic properties of soils using particle-size distribution and bulk density data. In: Indirect Methods for Estimating the Hydraulic Properties of Unsaturated Soils. University of California, Riverside, pp. 317–328. Churaev, N., Derjaguin, B., 1985. Inclusion of structural forces in the theory of stability of colloids and films. J. Colloid Interf. Sci. 103, 542–553. Conca, J., Wright, J., 1992. Diffusion and flow in gravel, soil, and whole rock. Hydrogeol. J. 1, 5–24. Derjaguin, B., Churaev, N., 1978. On the question of determining the concept of disjoining pressure and its role in the equilibrium and flow of thin films. J. Colloid Interf. Sci. 66, 389–398. Dormieux, L., Kondo, D., Ulm, F., 2006. Microporomechanics. Wiley, New York. Dullien, F.A., et al., 1992. Porous media: fluid transport and pore structure. Academic Press, San Diego. Dullien, F.A., Zarcone, C., Macdonald, I.F., Collins, A., Bochard, R.D., 1989. The effects of surface roughness on the capillary pressure curves and the heights of capillary rise in glass bead packs. J. Colloid Interf. Sci. 127, 362–372. Eshelby, J.D., 1957. The determination of the elastic field of an ellipsoidal inclusion, and related problems. In: Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. The Royal Society, pp. 376–396. Fredlund, D.G., Xing, A., 1994. Equations for the soil-water characteristic curve. Can. Geotech. J. 31, 521–532. Fredlund, M.D., Fredlund, D.G., Wilson, G., 1997. Prediction of the soil-water characteristic curve from grain-size distribution and volume-mass properties. In: Proceedings of the 3rd Brazilian Symposium on Unsaturated Soils, Rio de Janeiro, pp. 13– 23. Gonçalvès, J., Rousseau-Gueutin, P., de Marsily, G., Cosenza, P., Violette, S., 2010. What is the significance of pore pressure in a saturated shale layer? Water Resour. Res. 46, W04514. Grathwohl, P., 1998. Diffusion in Natural Porous Media: Contaminant Transport, Sorption/desorption and Dissolution Kinetics. Kluver Academic Publishers, Boston. Gruescu, C., Giraud, A., Homand, F., Kondo, D., Do, D., 2007. Effective thermal conductivity of partially saturated porous rocks. Int. J. Solids Struct. 44, 811–833. Hamamoto, S., Moldrup, P., Kawamoto, K., Komatsu, T., 2010. Excluded-volume expansion of archie’s law for gas and solute diffusivities and electrical and thermal conductivities in variably saturated porous media. Water Resour. Res. 46, W06514. Han, M., Youssef, S., Rosenberg, E., Fleury, M., Levitz, P., 2009. Deviation from archie law in partially saturated porous media: wetting film versus disconnectedness of the conducting phase. Phys. Rev. E 79, 031127. Hu, Q., Wang, J., 2003. Aqueous-phase diffusion in unsaturated geologic media: a review. Crit. Rev. Environ. Sci. Technol. 33, 275–297. Israelachvili, J., 1991. Intermolecular and Surface Forces. Academic Press, London. Kibbey, T.C., 2013. The configuration of water on rough natural surfaces: Implications for understanding air-water interfacial area, film thickness, and imaging resolution. Water Resour. Res. 49, 4765–4774. Kleinfelter-Domelle, N., Cushman, J., 2012. The role of connectivity in the theory of saturated/unsaturated flow through swelling media. Water Resour. Res. 48, W01543.

R. Yang et al. / International Journal of Multiphase Flow 79 (2016) 1–9 Koelemeijer, P.J., Peach, C.J., Spiers, C.J., 2012. Surface diffusivity of cleaved nacl crystals as a function of humidity: Impedance spectroscopy measurements and implications for crack healing in rock salt. J. Geophys. Res. 117, B01205. Koelman, J., De Kuijper, A., 1997. An effective medium model for the electric conductivity of an n-component anisotropic and percolating mixture. Physica A 247, 10– 22. Lemarchand, E., Davy, C.A., Dormieux, L., Chen, W., Skoczylas, F., 2009. Micromechanics contribution to coupled transport and mechanical properties of fractured geomaterials. Transport Porous Med. 79, 335–358. Lim, P., Barbour, S., Fredlund, D., 1998. The influence of degree of saturation on the coefficient of aqueous diffusion. Can. Geotech. J. 35, 811–827. Majumdar, A., Mezic, I., 1999. Instability of ultra-thin water films and the mechanism of droplet formation on hydrophilic surfaces. ASME Trans. J. Heat Transfer 121, 964– 971. Mehta, B.K., Shiozawa, S., Nakano, M., 1995. Measurement of molecular diffusion of salt in unsaturated soils. Soil Sci. 159, 115–121. Nakashima, S., 1995. Diffusivity of ions in pore water as a quantitative basis for rock deformation rate estimates. Tectonophysics 245, 185–203. Nguyen, S., 2014. Micromechanical approach for electrical resistivity and conductivity of sandstone. J. Appl. Geophys. 111, 135–140. Olesen, T., Moldrup, P., Gamst, J., 1999. Solute diffusion and adsorption in six soils along a soil texture gradient. Soil Sci. Soc. Am. J. 63, 519–524. Or, D., Tuller, M., et al., 2000. Flow in unsaturated fractured porous media: Hydraulic conductivity of rough surfaces. Water Resour. Res. 36, 1165–1177.

9

Raoof, A., Hassanizadeh, S., 2013. Saturation-dependent solute dispersivity in porous media: pore-scale processes. Water Resour. Res. 49, 1943–1951. Romkens, M., Bruce, R., 1964. Nitrate diffusivity in relation to moisture content of nonadsorbing porous media. Soil Sci. 98, 332–337. Sumner, A., Menke, E., Dubowski, Y., Newberg, J., Penner, R., Hemminger, J., Wingen, L., Brauers, T., Finlayson-Pitts, B., 2004. The nature of water on surfaces of laboratory systems and implications for heterogeneous chemistry in the troposphere. Phys. Chem. Chem. Phys. 6, 604–613. Tokunaga, T., 2011. Physicochemical controls on adsorbed water film thickness in unsaturated geological media. Water Resour. Res. 47, W08514. Wang, L., Bornert, M., Chanchole, S., Yang, D., Héripré, E., Tanguy, A., Caldemaison, D., 2013. Micro-scale experimental investigation of the swelling anisotropy of the callovo-oxfordian argillaceous rock. Clay Miner. 48, 391–402. Wang, L., Bornert, M., Héripré, E., Yang, D., Chanchole, S., 2014. Irreversible deformation and damage in argillaceous rocks induced by wetting/drying. J. Appl. Geophys. 107, 108–118. Weiss, J., Snyder, K., Bullard, J., Bentz, D., 2012. Using a saturation function to interpret the electrical properties of partially saturated concrete. J. Mater. Civil Eng. 25, 1097–1106. Yang, R.W., 2013. Contributions to micromechanical modelling of transport and freezing phenomena within unsaturated porous media (Ph.D. thesis). Université ParisEst.