Advances in Water Resources 97 (2016) 219–223
Contents lists available at ScienceDirect
Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres
Short communication
Prediction of capillary air-liquid interfacial area vs. saturation function from relationship between capillary pressure and water saturation E. Diamantopoulos a,∗, W. Durner b, T. Harter a a b
Department of Land, Air, and Water Resources, University of California, Davis, CA, USA Institute of Geoecology, TU Braunschweig, Braunschweig, Germany
a r t i c l e
i n f o
Article history: Received 6 May 2016 Revised 22 September 2016 Accepted 23 September 2016 Available online 23 September 2016 Keywords: Air-water interfacial area Saturation Contact angle Porous media Unsaturated flow
a b s t r a c t This short communication investigates if the capillary air-liquid interfacial area vs. saturation relationship Alv (S) can be predicted from the capillary pressure vs. saturation relationship S(h), using the theoretical sample scale model of Diamantopoulos and Durner (2013, 2015). We selected three published experimental datasets, where S(h) and Alv (S) relationships had been measured for the same porous media. The sample scale model was fitted to the retention curve S(h) of each porous medium and then used to predict the air-liquid interfacial area Alv (S). We also included in the analysis the thermodynamic models of Leverett (1941) and Grant and Gerhard (2007). For two sandy materials and especially for high saturation values, the model predicted the capillary Alv (S) successfully, which was in one case given by a porenetwork model simulation (Kibbey and Chen, 2012) and in the other case experimentally determined (Brusseau et al., 2006). For glass bead experiments, the contact angle needed to be fitted to properly describe the experimental Alv (S) curve. © 2016 Published by Elsevier Ltd.
1. Introduction The air-water interfacial area (Alv ) is an important parameter for characterizing phase distribution in unsaturated porous media. Hassanizadeh and Gray (1993), based on theoretical considerations, pointed out that different pore-scale air-water configurations may appear for the same porous medium at a given water content or saturation value. For this reason, an accurate description of hydraulic properties and multiphase flow in porous media must account for the geometry of the air-water interfaces (Chen et al., 2007a,b). This highlights the importance of the Alv vs. liquid saturation (S) relationship, which is ignored within the more commonly applied single-phase unsaturated flow models (Niessner et al., 2005). Determination of Alv is a non-trivial task. Generally, two types of methods are applied to determine Alv (S): tracer-based methods (Brusseau et al. 1997; Kim et al. 1997; Costanza-Robinson and Brusseau 2002) and imaging-based methods (Cheng et al., 2004; Culligan et al. 2004; Brusseau et al. 2006; Chen et al. 2007a; Porter et al. 2009). Both techniques have been widely used for determination of the Alv (S) relationship. However, recent discussions debate
∗
Corresponding author E-mail addresses:
[email protected],
[email protected] (E. Diamantopoulos). http://dx.doi.org/10.1016/j.advwatres.2016.09.012 0309-1708/© 2016 Published by Elsevier Ltd.
the type of interfaces that are measured by the two methods, respectively. It is generally agreed that tracer methods measure both, the capillary air-water interface and the film air-water interface, whereas the imaging techniques measure only capillary air-water interface area, Alv (Brusseau et al. 2006). An alternative method to investigate the relationship between interfacial area and saturation is the use of synthetic porous media (Karadimitriou and Hassanizadeh, 2012). A few studies exist, which investigate the Alv (S) relationship in synthetic micro-models (Chen et al., 2007; Pyrak-Nolte et al. 2008). Diamantopoulos and Durner (2013) developed a pore scale analysis for liquid retention and unsaturated hydraulic conductivity in angular capillary tubes as a function of the contact angle. Their analysis was based on the work of Tuller et al. (1999) and Or and Tuller (1999) without the consideration of films. For upscaling from the pore scale to the sample scale they assumed that the pore-size distribution can be represented by a special Gamma density function. Diamantopoulos and Durner (2015) increased the flexibility of their model by alternatively assuming a lognormal pore-size distribution. Their work provides analytical solutions for liquid retention, saturated/unsaturated hydraulic conductivity, liquid-air interfacial area, and specific surface area at the sample scale. The objective of this short communication is to test the model of Diamantopoulos and Durner (2013, 2015) by predicting the capillary air-water interfacial area vs. saturation relationship from experimental data of S(h) during primary drainage (h: pressure
220
E. Diamantopoulos et al. / Advances in Water Resources 97 (2016) 219–223
Table 1 Porosity, measured grain size and estimated retention curve parameters S(h) data for all materials using the lognormal model of Diamantopoulos and Durner (2015). The contact angle was assumed to be zero for Chen and Kibbey (2006) and Brusseau et al. (2006). For Culligan et al. (2004), the first case also assumes that the contact angle is zero, while the second case used Alv (S) data to fit the contact angle (see Fig. 3). Citation
Material
Porosity (ϕ )
Chen and Kibbey (2006) Kibbey and Chen (2012)
natural sand
0.362
Brusseau et al. (2006) Culligan et al. (2004)
natural sand glass beads
0.390 0.340
Weight percent grain size distribution (mm) 1%: 0.425–0.300 4%: 0.300–0.212 18%: 0.212–0.150 44%: 0.150–0.106 25%: 0.106–0.075 8%: 0.075–0.053 mean d: 0.234 mm 30%: 1–1.4 35%: 0.850 35%: 0.600
Lm (cm/-)
σ (-)
θ (°)
0.006
0.19
–
0.009 0.07 0.05
0.9 0.16 0.18
– – 47
head). Published data of Kibbey and Chen (2012), Brusseau et al. (2006) and Culligan et al. (2004) were used to predict Alv (S) to compare it with those measured in these works. For comparison, the analysis also applies the thermodynamic models of Leverett (1941) and Grant and Gerhard (2007).
on porous medium characteristics. Each of the triangular pores can be either fully saturated or partially saturated, dependent on the pressure head h. The macroscopic water saturation of a sample S(h) is calculated by summing up the water saturations in fully saturated and partially saturated pores:
2. Materials and methods
S ( h ) = S1 ( h ) + S2 ( h )
2.1. Experimental data
The equations for S1 and S2 with lognormal distributed poreside lengths (Eq. (1)) are:
Chen and Kibbey (2006) measured S(h) and Alv (S) relationships for a fine quartz sand (F-110, Chen et al., 2007b) by using five different concentrations of sodium octybenzene sulfonate (SOBS) as the interfacial tracer. In a later study, Kibbey and Chen (2012) showed that this method calculates both capillary and film (total) interfacial area. To do this, Kibbey and Chen (2012) calibrated a pore-network model (PNM) against experimental S(h) relationship for the same fine quartz sand. Afterwards, based on their model, they calculated total and capillary Alv (S) relationship. The former was compared against the experimentally measured Alv (S). Brusseau et al. (2006) measured S(h) and Alv (S) relationships for a natural sandy material by using two methods: gas phase partitioning tracer test and synchrotron X-ray microtomography. The first method measures both capillary and film (total) interfacial area whereas the second method can be used for the estimation of the capillary-associated Alv (S) relationship. Similarly, Culligan et al. (2004) used X-ray microtomography to estimate the S(h) and Alv (S) relationships for a glass bead column. We digitized the experimental data describing the primary drainage experiments by using the Web Plot Digitizer (Rohatgi 2012) on the published figures. Table 1 shows the experimental details for all three studies.
S1 ( h ) =
f (L ) =
Lσ
1 √
(ln( LmL ) )
2
2π
e−
2σ 2
, L > 0.
(1)
Parameter Lm is the median of the variable L, and σ is the standard deviation of the natural logarithm of L. Their values depend
L1 (h) L ln
L
√m
(3)
σ 2
1 −2 2σ 2 L e 1 + er f 2 m
S2 ( h ) = Bn
Lmin
L
2σ √ +√ σ 2 2
ln
Lmax
Lm
(4) L1 ( h )
where 2
Bn =
4r ( h ) Sn cot2
180 .
(5)
n
S1 (h) is the sample saturation of all fully saturated pores with a side length L greater than the minimum side length Lmin and smaller than L1 (h). S2 (h) is the sample saturation of the partially filled pores with L1 (h) < L < Lmax . Lmax stands for the pore sidelength of the biggest pore, which can be calculated from the soil’s air entry value or fixed to a relatively high value. The sample-scale liquid-vapor interfacial area is given by:
Al v ( h ) =
n(180 − 2(α + θ ) )π r (h ) 180An
1 −2 2σ 2 × L e 1 + er f 2 m
2.2. Diamantopoulos and Durner (2015) model The essentials of the model which describes both sample scale (as opposed to pore scale) S(h) and Alv (S) are presented by Diamantopoulos and Durner (2013). Briefly, the porous medium can be described as a bundle of angular tubes. In this study the tubes are assumed to have a shape of a equilateral triangular with a side length L. Diamantopoulos and Durner (2015) presented an upscaling scheme that assumed that the sample’s pore-size distribution, expressed by the probability distribution of the sample scale pore side-length, L, can be described by a log-normal distribution, f (L):
1 1 + er f 2
(2)
L
2σ √ +√ σ 2 2
ln
Lmax
Lm
(6) L1 ( h )
The variable definition for Eqs. (1)–(6) is given in Appendix A. Note that in our model only the interfaces between the corner fluid and the air phase contribute to the interfacial area, and that Alv (S) depends on the contact angle, θ , between the air-liquid interface and the solid surface at the pore wall. 2.3. Thermodynamic interfacial area model According to the thermodynamic approach, changes in interfacial area reflect the mechanical work done on the system, the magnitude of which is proportional to the area under the capillary pressure ((Pc ) [M T−2 L−1 ]) - saturation relationship (Grant and Gerhard, 2007). Leverett (1941) defined the total interfacial area per unit volume of porous media an [L2 L−3 ] as:
an = −ϕ
1
σwa
1
S
Pc (S )dS
(7)
E. Diamantopoulos et al. / Advances in Water Resources 97 (2016) 219–223
221
Fig. 1. a: Observed (grey circles) and fitted (red line) S vs. pF function for a sandy material (Chen et al., 2007b). pF is defined by pF = log10(h) with h in cm. b: Observed total Alv (S) (grey circles, Chen and Kibbey, 2012) and simulated capillary Alv (S) (grey triangles, Chen and Kibbey, 2012). Predicted an , awn and Alv (S) from Eqs. 9, 10 and 6 respectively.
where ϕ is the porosity, σ wa is the interfacial tension between water-air (this study) and the integral term depicts the area under the capillary pressure saturation curve. Eq. (7) is valid for the case of primary drainage under the assumption that S = 1 when Pc = 0. The capillary pressure is related to h by
Durner (2015) and the fitting was conducted by using the SCEUA algorithm (Duan et al. 1992). In a second step, we used Eq. (6) to predict Alv (S). Finally, we calculated both an and awn from Eqs. (9) and (10) respectively.
Pc = ρ gh
3. Results
(8)
where ρ is the density of water and g is the gravitational field [L T−2 ]. By substitution of Eq. (8) into (7) the latter becomes [M L−3 ]
an = −ϕ
ρg σwa
S 1
h ( S )d S
(9)
The total interfacial area an predicted by Eq. (9) represents the upper limit of interfacial area at any S value, and any model which predicts higher values of interfacial area is thermodynamically wrong/questionable. In this study we are interested in the capillary interfacial area and for this reason we calculated effective specific interfacial area awn based on the approach of Grant and Gerhard (2007). awn [L2 L−3 ] depicts the interfacial area that is only in contact between the bulk water phase and the air phase and is given by:
awn = −ψ (S ) Ed
ϕ
ρg σwa
1
S
h(S )dS
(10)
where ψ (S) [-] is the effective to the total specific/total specific interfacial area ratio and Ed is a term depicting energy dissipation effects and was taken equal to 0.21 (Grant and Gerhard, 2007). This value was adjusted based on the experimental data of Brusseau et al. (2006) for the case of primary drainage.The term ψ (Sw ) is approximated by (Grant and Gerhard, 2007) as:
ψ (S ) = −4.8871(S )4 + 13.18(S )3 − 11.901(S )2 +4.6119(S ) − 0.0073
(11)
(awn /an )
This equation was fitted to vs. saturation data from a pore-network modelling study (Dalla et al., 2002). 2.4. Approach Firstly, we fitted Eq. (2) to the S(h) experimental relationship. We assume that the residual saturation is equal to 0 for low values of the matric potential (i.e., high values of |h|). The fitted parameters were Lm and σ which both characterize the pore side-length distribution of the porous medium. The contact angle θ was assumed to be negligibly small for the model of Diamantopoulos and
Fig. 1a shows the measured S(h) relationship (grey points) for fine sand F-110 in the study of Chen et al. (2007b). Generally, the model (Eq. (2)) describes the experimental data quite well (red line). The fitted values of the parameters can be found in Table 1. Fig. 1b shows the experimentally determined Alv (S) as it was measured in Chen and Kibbey (2006) and later corrected in Kibbey and Chen (2012). The measured total air water interfacial area (mm−1 ) (grey circles) increases as water saturation decreases. Fig. 1b shows also the calculated from the PNM capillary Alv (S). The simulated Alv (S) increases as water saturation decreases, but it reaches a maximum around S≈0.5. Fig. 1b also shows the total interfacial area an and the effective specific interfacial area awn as predicted by Eqs. (9) and (10) respectively. an represents the upper limit of interfacial area according to the thermodynamic theory of Leverett (1941). As Kibbey and Chen (2012) pointed out, an is in accordance with the measured interfacial area, since the dynamic-interface tracer depletion method (used in Chen and Kibbey, 2006) measures both capillary and film interfacial area. Eq. (10) predicts that awn is lower than the one calculated from PNM capillary interfacial area, especially for 0.5 < S ≤ 1. Lastly, Fig. 1b also shows the Alv (S) relationship as predicted by Eq. (6). The model predicts the linear increase of Alv and matches the calculated Alv (S) relationship from Kibbey and Chen (2012) very well from full saturation to a low saturation of S ≈ 0.5. For high saturation values, this can be attributed to the pore-network model resampling a capillary bundle model. The model also predicts that Alv is largest at intermediate saturation, as it does not include the film component of Alv . Fig. 2a shows experimental and fitted S(h) relationships from the work of Brusseau et al. (2006). We compared the predicted Alv (S) with the measured relationship from the X-ray microtomography. Again, the fitted model describes the experimental reasonably well. In this short communication we focus on the capillary-associated Alv (S) relationship, since it covers most of the water saturation domain. The Alv (S) predicted from Eq. (6) is very close to the experimental one (Fig. 2b), at least for wet to intermediate conditions, 0.4 < S < 1. In the dryer part of the curve, the
222
E. Diamantopoulos et al. / Advances in Water Resources 97 (2016) 219–223
Fig. 2. a: Observed (grey circles) and fitted (red line) S vs. pF function for a natural sandy material (Brusseau et al., 2006). pF is defined by pF = log10(h) with h in cm. b: Observed (grey circles) and predicted an , awn and Alv (S) from Eqs. 9, 10 and 6 respectively.
Fig. 3. a: Observed (grey circles) and fitted S vs. pF function for a glass bead column (Culligan et al., 2004). pF is defined by pF = log10(h) with h in cm. b: Observed (grey circles) and predicted an , awn and Alv (S) from Eqs. 9, 10 and 6 respectively. The dashed red line assumes a contact angle of 0°, whereas the continuous red line is obtained by fitting Alv (S) by varying θ , yielding a fitted contact angle of 47°. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
model predicts a further linear increase in the Alv (S) curve. This deviation can be a result of the underlying model assumption (capillary bundle) but might be also due to limitations of the measurement method. The model of Grant and Gerhard (2007) gives a similar prediction to the model of Diamantopoulos and Durner (2015) but it requires the adjustment of parameter Ed to the experimental data. The experimental and measured S(h) relationships from the glass-bead study of Culligan et al. (2004) can be seen in Fig. 3a. The dashed red line shows the functions, if we fit of the retention curve and predict Alv (S), assuming that the contact angle θ is equal to 0 (Fig. 3b). Whereas S(h) is well described, the corresponding prediction clearly overestimated the measured Alv (S) function (Fig. 3b) by more than a factor of 2. As an alternative, we treated the contact angle θ in the function Alv (S) as additional free parameter in the inverse procedure and thus fitted both S(h) and Alv (S) simultaneously (continuous red line). With a contact angle of 47°, the functions describe the measured data very well. The required higher value of the contact angle for describing the experimental Alv (S) function can be attributed to the material, since for this study, it was glass beads which result to a finite value of the contact angle between water and the glass-beads (Lu et al. 1994). Garcia et al. (2016) used a simulated annealing method and conducted numerical simulations at the pore scale by minimizing the total interfacial energy in a hypothetical porous medium. Their
results also showed that the contact angle and porosity are the most important properties which control the extent of fluid interfaces. Also, the awn (S) predicted from the Grant and Gerhard (2007) model predicts very well the experimental data, without any need for further adjustment. 4. Conclusions In this short communication we investigated if Alv (S) can be predicted from S(h) at the sample scale. For three data sets from the literature with experimentally measured S(h) and Alv (S), we fitted the model of Diamantopoulos and Durner (2015) to S(h) and predicted Alv (S) for the same porous medium. We calculated also the total interfacial area an and the effective specific interfacial area awn from the work of Grant and Gerhard (2007). For the two sandy materials the model of Diamantopoulos and Durner (2015) was able to accurately predict the capillary Alv (S) from Kibbey and Chen (2012) that was simulated from PNM, and the experimental capillary Alv (S) from Brusseau et al. (2006), especially for high saturation values. In the case of glass beads, a nonzero value of the contact angle was required for the model of Diamantopoulos and Durner (2015) to also describe the experimental Alv (S) curve with sufficient accuracy. The model of Grant and Gerhard (2007) underestimated the measured interfacial area from Kibbey and Chen (2012) but described the measured interfa-
E. Diamantopoulos et al. / Advances in Water Resources 97 (2016) 219–223
cial area from Brusseau et al. (2006) very well. For the glass beads, it gave a very good prediction. To conclude, this study shows that both approaches of idealized pore space models and thermodynamic models may be used to predict the Alv (S) relationship from knowledge of the pore-size distribution, expressed by the S(h) relationship. Acknowledgments We thank Professor Dr. Tohren C. Kibbey for providing the experimental data used for model validation in Fig. 1a and b and for his constructive comments when reviewing this paper. Appendix A List of variables (Diamantopoulos and Durner, 2013; 2015) not described in the main text. In this study we assumed that the pores cross-section have the shape of equilateral triangle (n = 3, α = 30) Symbol
Dimensions
r(h)
L
Sn
– – – – –
α θ
n An
Definition Radius of the arc meniscus in the plane of the tube cross section θ cos (α + θ ) − π (1 − α +θ )] = tanα [ cos 2 90 sinα Half of the corner angle Contact angle Number of corners = n4 cot( πn )
References Brusseau, M.L., Popovicova, J., Silva, J.A.K., 1997. Characterizing gas-water interfacial and bulk-water partitioning for gas phase transport of organic contaminants in unsaturated porous media. Environ. Sci. Technol. 31, 1645–1649. Brusseau, M.L., Peng, S., Schnaar, G., Costanza-Robinson, M.S., 2006. Relationships among air-water interfacial area, capillary pressure, and water saturation for a sandy porous medium. Water Resour. Res. 42, W03501. http://dx.doi.org/10. 1029/20 05WR0 04058. Chen, D., Pyrak-Nolte, L.J., Griffin, J., Giordano, N.J., 2007a. Measurement of interfacial area per volume for drainage and imbibition. Water Resour. Res. 43, W12504. http://dx.doi.org/10.1029/20 07WR0 06021. Chen, L., Miller, G.A., Kibbey, T.C.G., 2007b. Rapid pseudo-static measurements of hysteretic capillary pressure-saturation relathionships in unconsolidated porous media. Geotech. Test. J. 30 (6), 1–10. Chen, L., Kibbey, T.C.G., 2006. Measurement of air–water interfacial area for multiple hysteretic drainage curves in an unsaturated fine sand. Langmuir ACS J. Surf. Colloids 22 (2006), 6874–6880. http://dx.doi.org/10.1021/la053521e. Cheng, J.-T., Pyrak-Nolte, L.J., Nolte, D.D., Giordano, N.J., 2004. Linking pressure and saturation through interfacial areas in porous media. Geophys. Res. Lett. 31, L08502. http://dx.doi.org/10.1029/2003GL019282.
223
Costanza-Robinson, M.S., Brusseau, M.L., 2002. Air–water interfacial areas in unsaturated soils: evaluation of interfacial domains. Water Resour. Res. 38, 1195. http://dx.doi.org/10.1029/20 01wr0 0 0738. Culligan, K.A., Wildenschild, D., Christensen, B.S.B., Gray, W., Rivers, M.L., Tompson, A.F.B., 2004. Interfacial area measurements for unsaturated flow through a porous medium. Water Resour. Res. 40, W12413. http://dx.doi.org/10.1029/ 20 04WR0 03278. Dalla, E., Hilpert, M., Miller, C.T., 2002. Computation of the interfacial area for two-fluid porous medium systems. J. Contam. Hydrol. 56, 25–48. Diamantopoulos, E., Durner, W., 2013. Physically-based model of soil hydraulic properties accounting for variable contact angle and its effect on hysteresis. Adv. Water Resour. 59, 169–180. http://dx.doi.org/10.1016/j.advwatres.2013.06.005. Diamantopoulos, E., Durner, W., 2015. Closed-form model for hydraulic properties based on angular pores with lognormal size distribution. Vadose Zone J. 14 (2). http://dx.doi.org/10.2136/vzj2014.07.0096. Duan, Q., Sorooshian, S., Gupta, V., 1992. Effective and efficient global optimization for conceptual rainfall-runoff models. Water Resour. Res. 28, 1015–1031. http:// dx.doi.org/10.1029/91WR02985. García, E.D., Boulet, P., Denoyel, R., Anquetil, J., Borda, G., Kuchta, B., 2016. Simulation of liquid–liquid interfaces in porous media. Colloids and Surf. A 496, 28–38. Grant, G.P., Gerhard, J.I., 2007. Simulating the dissolution of a complex dense nonaqueous phase liquid source zone: 1. Model to predict interfacial area. Water Resour. Res. 43, W12410. http://dx.doi.org/10.1029/20 07WR0 06038. Hassanissadeh, G.D., Gray, W.G., 1993. Thermodynamic basis of capillary pressure in porous media. Water Resour. Res. 29, 3389–3405. Karadimitriou, K.K, Hassanisadeh, S.M., 2012. A review of micromodels and their use in two-phase flow studies. Vadose Zone J. http://dx.doi.org/10.2136/vzj2011. 0072. Kibbey, T.C.G., Chen, L., 2012. A pore network model study of the fluid-fluid interfacial areas measured by dynamic-interface tracer depletion and miscible displacement water phase advective tracer methods. Water Resour. Res. 48, W10519. Kim, H., Rao, P.S.C., Annable, M.D., 1997. Determination of effective air-water interfacial area in partially saturated porous media using surfactant adsorption. Water Resour. Res. 33, 2705–2711. Leverett, M.C., 1941. Capillary behaviour in porous solids. Trans. AIME 142, 152–170. Lu, T.X., Biggar, J.W, Nielsen, D.R., 1994. Water movement in glass bead porous media. 1. Experiments of capillary rise and hysteresis. Water Resour. Res. 30 (12), 3275–3281. Niessner, J., Helmig, R., Jakobs, H., Roberts, J.E., 2005. Interface condition and linearization schemes in the Newton iterations for two-phase flow in heterogeneous porous media. Adv. Water Resour. 28, 671–687. Or, D, Tuller, M, 1999. Liquid retention and interfacial area in variably saturated porous media: upscaling from single-pore to sample-scale model. Water Resour. Res. 35 (12), 3591–3605. http://dx.doi.org/10.1029/1999WR900262. Porter, M.L., Schaap, M.G., Wildenschild, D., 2009. Lattice-Boltzmann simulations of the capillary pressure-saturation-interfacial area relationship for porous media. Adv. Water Resour. 32, 1632–1640. Pyrak-Nolte, L.J., Nolte, D.D., Chen, D., Giordano, N.J., 2008. Relating capillary pressure to interfacial areas. Water Resour. Res. 44, W06408. http://dx.doi.org/10. 1029/20 07WR0 06434. Rohatgi A., 2012. WebPlotDigitalizer: HTML5 based Based online Online tool Tool to extract Extract numerical Numerical data Data from plot Plot imagesImages. Version 3.9. URL http://arohatgi.info/WebPlotDigitizer/app/. Tuller, M., Or, D., Dudley, LM., 1999. Adsorption and capillary condensation in porous media: liquid retention and interfacial configurations in angular pores. Water Resour. Res. 35 (7), 1949–1964. http://dx.doi.org/10.1029/1999WR90 0 098.