Journal of Contaminant Hydrology 247 (2022) 103984

Contents lists available at ScienceDirect

Journal of Contaminant Hydrology
journal homepage: www.elsevier.com/locate/jconhyd

Comparison of methods to estimate air-water interfacial areas forevaluating PFAS transport in the vadose zone
Jeff A.K. Silva a, *, Jiří Šimůnek b, John E. McCray caArclight Research & Consulting, Golden, CO, USA
Department of Environmental Sciences, University of California-Riverside, USAc
Hydrologic Science and Engineering Program, Civil & Environmental Engineering Department, Colorado School of Mines, Golden, USAb

A R T I C L E I N F O

A B S T R A C T

Keywords:
PFAS
Air-water interfacial adsorption
Air-water interfacial area
Vadose zone transport
Soil texture and surface roughness effects
Modeling

When performing calculations or numerical simulations for the fate and transport of PFAS and other surfaceactive solutes in the vadose zone, accurately representing the relationship between the area of the air-waterinterfaces (Aaw) as a function of water saturation (Sw), and changes in that relationship resulting fromchanges in soil texture, are equally important as accurately characterizing interfacial adsorption coefficients andthe concentration dependence for PFAS solutes. This is true because the magnitude of the Aaw directly governs
the degree of air-water interfacial adsorption, which contributes to the transport retardation of these solutes
within unsaturated porous media. Herein, a well-known thermodynamic-based model for predicting the Aaw-Sw
relationship is evaluated through comparisons to literature data collected using various measurement techniques
for model sands and a limited number of soils using data collected from the current published literature. Thispredictive model, herein termed the Leverett thermodynamic model (LTM), relies on the characterization of thesoil-water retention curve (SWRC) for a given soil, using the van Genuchten (VG) equation for the pressure headvs-Sw relationship. Therefore, methods to estimate the VG equation parameters are also compared as to the AawSw relationships predicted. Comparisons suggest that the LTM provides the best estimate of the actual Aaw-Sw
relationships for water containing non-surface-active solutes. Because PFAS solutes are also surface-active, Aaw
measurement methods utilizing surface-active tracers are considered to provide the most accurate representationof the Aaw-Sw relationship for these solutes. Differences between Aaw-Sw relationships derived from tracer
methods and the LTM are described in relation to media surface roughness effects. Based on the available
literature data, a practical empirical model is proposed to adjust the LTM prediction to account for the effects ofsurface roughness on the magnitude of the Aaw for surface-active solutes. Finally, example retention calculations
are performed to demonstrate the sensitivity of the predicted Aaw-Sw relationship on the vadose zone transport ofof a representative PFAS, perfluorooctane sulfonate.

1. Introduction
Per- and polyfluoroalkyl substances (PFAS) are a class of emerging
contaminants that have received significant attention in recent years
due to human health concerns, their widespread use and distribution in
the environment, and their general resistance to degradation in theenvironment. PFAS groundwater contamination is often associated withthe leaching of these compounds from commercial products initiallydeposited in unlined landfills or on the ground surface (e.g., from landapplied bio-solids and past releases of aqueous film-forming foam
(AFFF) products at fire training areas). In these release scenarios, PFAS

must first transport through the vadose zone prior to reaching under­
lying aquifers. PFAS has been shown to accumulate within the shallow
vadose zone, serving as a long-term source of groundwater contamina­tion (Shin et al., 2011; Xiao et al., 2015; Weber et al., 2017; Anderson
et al., 2006; Anderson et al., 2019). Therefore, it is important to un­derstand and characterize PFAS fate and transport mechanisms within
the vadose zone to better assist environmental risk assessment andfuture remediation decision-making.
PFAS are surface-active chemicals (i.e., surfactants). As surfactants,
PFAS will accumulate at fluid interfaces during transport through thevadose zone. Of particular interest in this paper is the adsorption of PFAS* Corresponding author.
E-mail address: jsilva@arclightrc.com (J.A.K. Silva).
https://doi.org/10.1016/j.jconhyd.2022.103984
Received 19 July 2021; Received in revised form 22 February 2022; Accepted 27 February 2022
Available online 7 March 2022
0169-7722/© 2022 Elsevier B.V. All rights reserved.

J.A.K. Silva et al.

Journal of Contaminant Hydrology 247 (2022) 103984

at air-water interfaces (AWI), although other important interfacial
processes have been identified as additional contributors to overall PFAS
retention (e.g., Guelfo and Higgins, 2013; McKenzie et al., 2016; Brus­seau, 2018). PFAS adsorption at AWIs is significant under certain con­ditions (e.g., Brusseau, 2018; Silva et al., 2019; Guo et al., 2020; Silva
et al., 2020), depending on the specific surface activity of a given PFAS
and the magnitude of the area of the AWI (Aaw). Aaw is inversely related
to volumetric water content (θw) or water saturation (Sw = θw/porosity).
Variability in Aaw depends on the particle size distribution in watercontent, which in turn depends on textural properties of the vadose zone
soil. Each of these dependencies can vary spatially and temporally
within the vadose-zone environment, complicating transport
calculations.
Methods to characterize the Aaw-Sw relationship have included
theoretical models (based on empirical relationships or assumptions ofpore geometry [e.g., Leverett, 1941; Bradford and Leij, 1997; Gvirtzman
and Roberts, 1991; Or and Tuller, 1999; Oostrom et al., 2001; Jiang
et al., 2019]) and experimental measurement methods utilizing surfaceactive gas-phase and aqueous-phase chemical tracers (e.g., Kakare andFort, 1996; Kim et al., 1997; Silva, 1997; Kim et al., 1999; Schaefer et al.,
2000; Costanza-Robinson, 2001; Silva et al., 2002; Araujo et al., 2015;
Brusseau et al., 2020). More recently, direct imaging via x-ray micro­
tomography (XMT) methods have additionally been employed (e.g.,
Brusseau et al., 2007; Costanza-Robinson et al., 2008; Araujo andBrusseau, 2020). The various methods used for Aaw measurement
methods have been shown to provide different Aaw values for a given Sw
value, as well as different Aaw-Sw relationships, particularly at lower
moisture conditions. For example, XMT methods provide linear Aaw-Sw
relationships throughout the entire range of Sw. Gas-phase tracer
methods (i.e., those in which a gaseous or vaporous tracer is introduced
to a mobile air phase and where the aqueous-phase and Aaw are assumed
static) typically provide highly non-linear Aaw-Sw relationships and Aaw
values are considerably greater than those provided by other measure­ment methods.
With respect to aqueous tracer methods, there are two main types.
The first type involves introducing the interfacial tracer with aqueous
advection wherein the tracer contacts the available interfacial area
under steady unsaturated flow (e.g., Brusseau et al., 2015 and references
therein). For the second method, the tracer solution is introduced to an
unsaturated sand column and the Aaw is estimated via a mass balancecalculation from the equilibrium mass of tracer lost from the test column
when allowed to drain to a steady-state moisture condition (Schaefer
et al., 2000; Araujo et al., 2015). Schaefer et al. (2000) observed linear
Aaw-Sw relationships while Araujo et al. (2015), using similar methods,
observed a non-linear increasing trend in Aaw with increasing Sw. Bothgaseous and aqueous tracer methods provide Aaw values that are greaterthan those derived from XMT measurements at the same Sw condition
and the differences increase with decreasing Sw. The differences in Aaw
values provided by the different measurement methods has been argued
to result from each method characterizing different fractions of the total
Aaw (e.g., Aaw associated with capillary-held pendular water and Aaw
associated with thin water films adsorbed to mineral surfaces) (e.g.,
Costanza-Robinson, 2001; Brusseau et al., 2006).
The variability in Aaw-Sw relationships presented in the literature canlead to confusion over which relationships are most appropriate when
performing fate and transport calculations for PFAS or other surfaceactive solutes. The most commonly held view based on the collective
literature to date is to use Aaw-Sw relationships that best represent thetransport conditions being investigated. For PFAS of current environ­mental concern, transport within the vadose zone occurs in the aqueous
phase and Aaw-Sw relationships derived from aqueous-based methods
would therefore be most representative (Silva et al., 2020; Guo et al.,2020: Brusseau et al., 2020; Brusseau et al., 2019), as these methods are
considered to characterize Aaw accessible to surface-active solutes pre­sent in the aqueous phase. However, discrepancies between thesemethods remain, and too little data representative of a variety of soil

types are available to make a clear judgment on the utility of the Aaw-Sw
relationships provided. Likewise, the diversity of theoretical models
available in the literature to predict the Aaw-Sw relationship for a given
porous media adds to the complexity and uncertainty.
Herein, the efficacy of one of several theoretical models to predict
Aaw-Sw relationships, a modification of the thermodynamic-based model
first conceptualized by Leverett (1941), hereafter referred to as the
Leverett Thermodynamic Model (LTM), is evaluated. Data used in thisevaluation represents a variety of model porous media and natural sandy
soils, with the goal of improving the accuracy of transport calculations
for PFAS and other similarly surface-active solutes. This theoreticalmodel has been used to simulate colloid and bacteria transport (Šimůnek
et al., 2016) and was recently implemented in an additional modified
version of the HYDRUS unsaturated flow and transport model to simu­
late PFAS transport in the vadose zone (Silva et al., 2020). Predicted
results are compared to Aaw-Sw relationships derived from the literaturethat utilized XMT and aqueous tracer-based measurements. Aaw-Sw re­lationships predicted by the LTM method are directly dependent on thecharacteristics of the soil-water retention curve (SWRC) for a given
porous media. Here, the well-known van Genuchten (VG) equation was
selected to model SWRCs. Therefore, evaluation and discussion of the
impact of key VG equation parameters on predicted Aaw-Sw relationships
are additionally provided. Finally, an empirical relationship is proposedthat can be used to scale the LTM to be more representative of Aaw values
measured by surface-active solutes, thereby improving estimates of
PFAS retention within the vadose zone due to AWI adsorption.
2. Materials and methods
2.1. AWI area calculation
The method used here to estimate Aaw as it changes with soil watercontent (θ) was proposed by Leverett (1941) and implemented by
Bradford et al. (2015) as:∫∫
1 θs
ρ g θsAaw (θ) =Pc (θ)dθ = w
h(θ)dθ
(1a)
γo θ
γo θ
dAaw (Θ) Pc ρw gh
=
(1b)
=
γo
γo
dθ
where θs is the saturated water content (L3/L3), Pc is the capillarypressure (M/LT2), h is the pressure head (L), γ 0 is the air-water surface
tension (M/T2), ρw is the density of water (M/L3), and g is the gravita­
tional acceleration (L/T2). The integration provides Aaw in units of L2/L3(or L− 1), where the L3 term in the numerator is the total volume of air,water, and solid phases. The magnitude of Aaw increases with decreasing
θ and approaches the geometric surface area (GSA) of the porous media
as θ approaches the residual (irreducible) water content (θr).
This approach relies on a specific model characterizing the change inwater content in response to a change in h. Here, we selected the wellknown van Genuchten (VG) model (van Genuchten, 1980):⎧⎪⎨ θr +
θ(h) =⎪⎩

θs − θr[1 + |αh|n ]

/ h≤0

1− 1

n

(2)

θs h > 0
where θr has units of L3/L3), α (1/L) is a parameter related to the inverseof the air-entry pressure head for the porous medium, and n (− ) is aparameter related to the pore-size distribution. Eq. (1) describes the soilwater retention curve (SWRC) for a porous medium.
The use of Eq. (1) and similar methods to estimating Aaw-Sw re­lationships have been previously referred to as the “thermodynamic”
approach (e.g., Dalla et al., 2002; Grant and Gerhard, 2007; Schroth
et al., 2008;) because, as proposed by Leverett (1941), it relates changes
in Aaw to the thermodynamic work done within a unit volume of2

J.A.K. Silva et al.

Journal of Contaminant Hydrology 247 (2022) 103984

unsaturated porous media when changes in water contents occur due tochanges in capillary pressure. The integration of Eq. (1), i.e., the area
under the SWRC, yields the Aaw-Sw relationship. The endpoint of the Aaw
prediction is the θr for a given porous medium, where it is assumed thatsoil moisture exists as pendular water and as a thin adsorbed layer ofwater coating mineral surfaces, and advective flow does not occur. Thismodel is hereafter referred to as the Leverett thermodynamic model(LTM). The LTM is applicable strictly at the continuum (or representa­
tive elementary volume (REV)) scale, similar to the SWRC. Therefore,
media characteristics not captured by the SWRC would not be repre­
sented in the Aaw predicted by the LTM.

2.4. Estimating VG model α and n parameters2.2. X-ray microtomography dataset

where α1d is the α for Cu = 1, Nα is α normalized by α1d, and α corrected
for Cu is α1dNα. d60 is the particle diameter (L) of which 60% of theparticles in the soil sample are finer. Similarly, n (drying curve) is
calculated as:

VG model parameters α and n for each media type evaluated in thiswork were developed using the pedotransfer function presented by
Benson et al. (2014) (hereafter referred to as the Benson PTF), whichallows correction of VG parameters α and n for changes in the particle
size distribution based on changes in the uniformity coefficient (Cu). For
α as 1/kPa (drying curve case, from Benson et al., 2014):

The results of the LTM predictions are compared to those presented
by Costanza-Robinson et al. (2008), who reported Aaw-Sw relationships
developed from synchrotron x-ray microtomography (XMT) measure­ments for glass bead media and model sands. The XMT method is an
imaging methodology that allows for the direct calculation of the Aaw-Sw
relationships. The current interpretation of the XMT method results isthat the Aawvalues provided represent the combined contributions of theareas associated with capillary water (e.g., pendular rings at points ofcontact between mineral particles) and adsorbed water films on wetted
surfaces within the pores (e.g., Brusseau et al., 2007; Costanza-Robinson
et al., 2008). XMT cannot resolve interfacial area associated withmicroscale surface heterogeneity (e.g., surface roughness), when pre­sent, and thus provides maximal Aaw values that are nearly equivalent tothe calculated smooth specific surface area of the porous media (Brus­seau et al., 2007; Costanza-Robinson et al., 2008).
As presented by Costanza-Robinson et al. (2008), XMT provided
linear Aaw-Sw relationship models for each media type that had thefollowing form (after Araujo and Brusseau, 2020):
Aaw = GSA(1 − Sw )(5)Nα = 0.99Cu − 0.54(6)n1d = 14exp( − 0.434d60 )(7)

Nn (Cu > 2.2) = − 0.0033Cu + 0.379(8)where n1d is the n for Cu = 1, Nn is n normalized by n1d, and n corrected
for Cu is n1dNα. The Benson PTF utilizes regressions based on d60 for the
media (i.e., the soil particle diameter for which 60% of the total are finer
by weight). The d60 values for the media used here were not reported;
however, the d50 values were. For reasonably uniform porous media, thenear equivalence of d60 and d50 values can be and was assumed herein.
3. Results and discussion
3.1. Porous media properties
Porous media properties and VG parameters compiled and calculated
for this evaluation are provided in Table 1. For the various glass bead
media, calculated VG parameters α and n are consistent with bothmeasured values as well as values calculated using alternative theoret­
ical packing arrangements of spherical particles (Benson et al., 2014;
Sweijen et al., 2017). Likewise, calculated α and n are consistent with
reported measured values for the Accusand media (Schroth et al., 1996).
According to the Benson PTF, α values follow a linearly increasing trendwith increasing values of the d60 for these media. Calculated n values
trend as a decreasing power-law function with increasing Cu. For theVinton soil, α an n values show some variability for the different
parameter estimation methods used (i.e., coefficients of variation (CV)
= 25% and 15% for the estimates of α and n, respectively). The esti­
mated values are reasonable for this type of soil.

(3)

where GSA is the geometric surface area of the medium (L2/L3) based on
the smoothed-sphere assumption (i.e., the medium is composed of acollection of smooth spheres of an identical diameter), calculated as:
GSA = 6(1 − n)/d50

α1d = 1.34d60(4)where n is the media porosity (L3/L3) and d50 is the median particle
diameter (L). Costanza-Robinson et al. (2008) actually derived thefollowing equation: Aaw = SA(0.9031 − 0.9112Sw). However, throughour analysis, we have found Eq. (4) to be a reasonable simplification torepresent the measured data presented by Costanza-Robinson et al.
(2008). Hereafter, Eq. (3) will be used to represent the XMT-measured
data for the soils used in this work.

3.2. VG parameters α and n and the effect on Aaw prediction

2.3. Additional soil properties data

The effect of soil texture on the values of the VG equation α and nparameters used to represent the SWRC is well known. For clean sands
and other coarse granular porous media, there is a strong relationshipbetween particle/pore-size distribution and the character of the SWRC
(Yang et al., 1994; Benson et al., 2014). For uniform porous media, adecrease in α relates to an increase in the air entry pressure/suction,with a slight increase in the slope of the SWRC (i.e., a decrease in themagnitude of n) between θs and θr as the particle size becomes finer. As
the particle-size distribution becomes increasingly non-uniform, airentry pressures, and thus the value of α, will vary depending on thenature of the particle-size distribution and corresponding changes inparticle packing arrangements. However, it is the slope of the SWRC,
and thus the magnitude of n, that is most affected by the degree ofparticle-size non-uniformity.
When using Eq. (1) and Eq. (2) to describe the Aaw-Sw relationship,
the impacts of α and n on the shape of the SWRC are described below.
SWRCs and their corresponding Aaw-Sw relationships calculated for the

Additional soil properties data were derived from literature sources
for two natural sands (Torrens et al., 1998; Costanza-Robinson, 2001;
Brusseau et al., 2006; Peng and Brusseau, 2012) and used for compari­son of LTM predictions. For example, measured VG parameters for twodesert soils (i.e., Vinton sandy soil and Hayhook sandy loam) were ob­tained from the literature (Peng and Brusseau, 2012) and directly used
to develop Aaw-Sw relationships via the LTM approach. Soil textural data
(i.e., percentages of sand, silt, and clay) were also derived from theliterature for the same Vinton soil (Costanza-Robinson, 2001) andHayhook soil (Torrens et al., 1998) that were used to develop VG pa­
rameters via the well-known Rosetta Lite (v. 1.1) neural network modelusing the SSCBD (or H3) model described by Schaap et al. (2001), which
included the bulk density (ρb) of the media as an additional predictive
parameter. Finally, the LTM approach was applied to VG parameter data
(i.e., specifically α and n values) derived from the pedotransfer function(PTF) presented by Benson et al. (2014), as described below.
3

J.A.K. Silva et al.

Journal of Contaminant Hydrology 247 (2022) 103984

Table 1
Porous media properties and VG parameters collected or calculated for this evaluation.
Porous media characteristics
Reference for porous media characteristics

Costanza-Robinson et al. (2008)

Costanza-Robinson et al. (2008); Peng andBrusseau, 2005, 2012)
Araujo and Brusseau (2020)

Media
Glass Beads
Fine
MediumCoarseMixtureAccusand Sands
MediumCoarseGranusil Sands
Fine
CoarseMixtureVinton Sandy Soil
Vinton (Benson
PTF)
Vinton(measured)dVinton (Rosetta)e
Hayhook Loamy
Sand
Hayhook (Benson
PTF)
Hayhook(measured)dHayhook(Rosetta)e
1.1 mm Glass
Beads

a

VG equation parameters

d50
(mm)

d10
(mm)

Cu

GSA
(cm− 1)SSABETb
− 1

0.099
0.345
0.650
0.403

0.090
0.288
0.591
0.088

1.1
1.2
1.1
4.6

365
113
55
97

780
211
296
545

0.375
0.516

0.313
0.430

1.2
1.2

101
73

0.200
0.513
0.211

0.118
0.302
0.073

1.7
1.7
2.9

–

–

0.234
–

θr

θs

α

n

2
2
5
6

0.05
0.02
0.02
0.02

0.398
0.350
0.404
0.348

0.012
0.041
0.081
0.023

13.0
11.1
10.3
4.40

1776
2103

18
29

0.02
0.02

0.369
0.372

0.045
0.061

10.91
10.26

162
65
163

8473
5648
5175

52
87
32

0.05
0.05
0.05

0.425
0.430
0.427

0.024
0.061
0.022

8.19
7.15
4.85

–

–

–

–

0.05

0.419

0.019

4.83

0.098
–

2.4
–

155
–

51,684
–

334
–

0.06
0.05

0.390
0.360

0.026
0.032

3.66
3.85

–

–

–

–

–

–

0.05

0.360

0.008

4.20

0.260

0.016

16

148

144,320

978

0.05

0.390

0.024

2.26

–

–

–

–

–

–

0.05

0.352

0.033

1.96

1.1

1.1

1

32

28

0.9

0.01

0.340

0.152

8.46

(cm

)

SRF

c

(cm− 1)

Note: Porous media characteristics presented were derived from the sources cited in the table (leftmost column). VG equation parameters were calculated using the
Benton PTF or as otherwise derived from sources identified in the table by footnote. aCu = d60/d10 or d50/d10 as an estimate in this case. breported specific surface areafrom N2/BET, where BET is the Brunauer-Emmett-Teller equation. cSRF is the shape/roughness factor = SSABET/GSA. dmeasured VG parameters as reported by
Brusseau et al. (2006). eVG parameters estimated using Rosetta program (Schaap et al., 2001) and the texture-bulk density model (97% sand, 1.8%silt, 1.2% clay, bulk
density = 1.63). ϴr values derived from Benson et al. (2014), Schroth et al. (1996), Sakaki and Illangasekare (2007), Brusseau et al. (2006). VG parameters for 1.1 mm
glass beads derived using Benson PTF.

glass beads and model sands used in this evaluation are presented inFig. 1. As demonstrated for the uniform and spherical glass beads(Fig. 1a and d), a decrease in particle size, and a corresponding decreasein α, results in an increase in the air-entry pressure for the SWRC and acorresponding increase in the overall slope of the Aaw-Sw relationship.
Mathematically, increased Aaw with decreasing particle size occurs dueto the overall increase in the integrated area under the SWRC. Mecha­
nistically, however, a reduction in particle size relates to an increase in
the GSA, and thus an overall increase in calculated Aaw for Sw < 1. As
demonstrated for the glass bead mixture, the reduction in n increases thedegree of curvature of the Aaw-Sw relationship as Sw decreases, which isalso the case for the corresponding SWRC.
The Aaw-Sw relationships for the Granusil sand mixture also exhibit
an increasing degree of curvature, albeit less than those observed for the
spherical glass beads. In this case, the Aaw values predicted are overall
greater than those predicted for the glass bead mixture despite the factthat the Cu for the glass bead mixture is a factor of 1.6 greater than thatof the Granusil mixture. The differences in these Aaw-Sw relationships are
the result of the differences in the d50 reported for these sand mixtures,
which impact the α and n values provided by the Benson PTF. From thisanalysis, it is observed that for uniform sandy porous media, the particle(pore) size – and thus the magnitude of α – plays a dominant role in themagnitude of Aaw values predicted by the LTM. As the particle size
distribution becomes increasingly less uniform (i.e., a wider particle size
distribution), the contribution of the n parameter to the predicted Aaw
values become increasingly important.

Costanza-Robinson et al. (2008), are represented in Fig. 1 to provide acomparison to the LTM predictions. For the glass beads, the Aaw-Sw re­lationships predicted by the LTM (using the α and n values provided by
the Benson PTF) are very similar to those provided by the XMT mea­
surements. Both methods practically predict the same maximum inter­
facial areas (i.e., the endpoint Aaw prediction at Sw = 0), which are
estimated using the smooth surface GSA values represented in Table 1.
For the fine glass bead media, the LTM prediction deviates upward fromthe XMT results with increasing separation as Sw approaches Sr, which
results from the lower value of α provided by the Benson PTF in responseto the decrease in particle size. Costanza-Robinson et al. (2008) describe
the XMT results as being less accurate for the fine bead media due to thepoor image segmentation quality that resulted in only a qualitative
estimation of Aaw values. However, both methods provided reasonably
similar Aaw-Sw relationships. This is true for the highly uniform andlargely spherical Accusand media as well (e.g., particle roundness = 0.8,
sphericity = 0.8–0.9 (Schroth et al., 1996; Hamamoto et al., 2016). The
observed similarity in Aaw values implies that both methods are char­acterizing the same effective interfacial areas for these uniform, spher­
ical, and smooth porous media.
For the glass bead mixture, the Aaw values predicted by the LTM areoverall greater than those by the XMT method. In this case, it appears
that the XMT may be underestimating Aaw values for this mixture.
Costanza-Robinson et al. (2008) noted that pore-water within the XMTimaged bead mixture sample was observed to be non-uniformly
distributed within the images analyzed due to preferential capillary
filling of small pores associated with the fine bead fraction (i.e., a pre­sented imaged section for the bead mixture sample used in the analysis
shows an area of water associated with the fine bead fraction that en­compasses roughly a third of the image area); whereas pore-water was

3.3. Comparing literature-derived XMT results to LTM predictions
The results of Eq. (3), based on the XMT results presented by
4

J.A.K. Silva et al.

Journal of Contaminant Hydrology 247 (2022) 103984

Fig. 1. SWRCs (a-c) and associated Aaw-Sw relationships (d-f) for porous media used in this evaluation. Linear dashed Aaw-Sw relationships are those estimated andreported by Costanza-Robinson et al. (2008) using XMT methods. Solid lines are the LTM predictions. Filled symbols are the calculated GSA values for the media.
Open symbols are measured data from Sakaki and Illangasekare (2007).

observed to be homogeneously distributed in the uniform bead media
sample images. The presence of heterogeneously distributed water in the
non-uniform sample is not surprising in and of itself. However, it calls
into question whether the Aaw value estimated from this sampled image
would be reproducible for another sampling at the same moisture con­dition and at the scale of measurement. This issue also raises the ques­tion as to whether the sample scale for the XMT measurement meets therequirement of a REV for this non-uniformly sized glass bead mixture.
Complications related to the REVs and XMT Aaw calculations for nonuniform media have been recently documented (Costanza-Robinson
et al., 2011; Wang et al., 2019). The increased slope of the LTMpredicted Aaw-Sw relationship is the result of the reduction in the αparameter provided for this non-uniform bead mixture (Cu = 4.6). For
the VG equation, parameters vary systematically with the median par­ticle size; α should decrease as the median particle size decreases and as

the breadth of the particle size distribution increases, and n should
decrease as the breadth of the particle size distribution increases (Ben­son et al., 2014).
Aaw values for the Granusil sands predicted by the LTM, again using α
and n values provided by the Benson PTF, were also overall larger thanthose for the XMT results for both the fine- and coarse-grained fractions
and the mixture. However, the reason for the difference in Aaw values, inthis case, does not appear to be related to the heterogeneous pore-water
distribution issues observed for the glass bead mixture. The pore-water
distribution in the XMT images presented by Costanza-Robinson et al.
(2008) appears to be more uniform within the Granusil sands in allcases. Here, we interpret the differences between XMT and LTM esti­mates as being related to the Granusil sands' angularity (roundness =0.2; sphericity = 0.7) (Hamamoto et al., 2016) and increased roughnessof the particle surfaces. As noted by Costanza-Robinson et al. (2008), the5

J.A.K. Silva et al.

Journal of Contaminant Hydrology 247 (2022) 103984

XMT method uses surface smoothing algorithms to minimize imaging
artifacts as a part of estimating Aaw values, which removes micro-scalesurface topographic heterogeneity (i.e., surface roughness) of media
particles. XMT, therefore, provides smoothed-surface estimates of Aaw
that compare well with GSA values. Conversely, previous research has
shown that particle surface roughness is captured in the results of media
fluid retention curves (Dullien et al., 1989; Zheng et al., 2015). With
respect to the SWRC, the angularity of the Granusil media allows greaterpacking density than for more rounded sands of the same particle size,
which results in an smaller pores and subsequently an increased airentry pressure for this media during drainage, with a corresponding
reduction in the VG α parameter. Likewise, an increased media surface
roughness would also manifest in a reduced α parameter. The net effectof these media characteristics would then be an increased slope of the
Aaw-Sw relationship predicted by the LTM than would be measured by
XMT, as demonstrated in Fig. 1f. Given the agreement between
measured data and the predicted SWRCs shown in Fig. 1c, the Benton
PTF appears to be providing VG parameters that are indirectly, but
appropriately, accounting for the effects of the angularity and surface
roughness of the Granusil sands.
3.4. Comparing XMT and LTM predictions for natural Sandy soils
The majority of the existing literature data relating to estimating Aaw
for unsaturated porous media is limited to some model sands (i.e, con­structed using well-sorted sand fractions) and a few sandy soils (i.e., theVinton and Hayhook soils described previously). Relevant data collected
for these sands are also presented in Table 1. To our knowledge, when
this paper was prepared, XMT measurement data had been only reported
for the Vinton soil (Brusseau et al., 2006). The LTM was applied using
VG parameters derived from the Benson PTF, published measured pa­
rameters, and those derived from the Rosetta model as described in
Section 2.2. For the Rosetta model, the textural classification of theVinton soil was 97% sand, 1.8% Silt, 1.2% clay, and ρd= 1.46 (CostanzaRobinson, 2001; Peng and Brusseau, 2005). Likewise, the Hayhook soil
was classified as containing 86% sand, 4.7% silt, and 8.8% clay (Torrens
et al., 1998), with ρb= 1.64 g cm− 3 (Peng and Brusseau, 2005).
The Aaw-Sw relationships developed for the Vinton and Hayhook soils(Fig. 2) demonstrate the dependence of the LTM prediction on the se­lection of α and n provided by different parameter estimation methods
selected. Note that for both soils, the Benson PTF predicts lower α values(Table 1), which translates to a higher air-entry pressure for the SWRCs
and increased Aaw values than those predicted when using the othermethods at the same Sw (Fig. 2). Additionally, the Benson PTF provides
values of n that are larger than those derived using the other methods. In
comparing the SWRCs derived from the Benson PTF and the measured
VG parameters for these natural sandy soils, the Benson PTF appears tooverestimate the moisture retention contribution of the fines fraction at
higher Sw. As noted by Benson et al. (2014), the PTF was developed
using sands for which the fine fraction (i.e., particles <0.075 mm) wasremoved. Therefore, exclusion of the fines fraction appears to beweighting the effect of Cu more heavily on the predicted value of α thanon the value of n when fines exist in the sample and for Cu > 2.
To correct this, the following simple modification to the Benson PTF
was made:Nα = 0.99Cu − 0.2 , and(9)Nn (Cu > 2.2) = − 0.0033Cu + 0.279(10)Fig. 2. SWRCs (a-b) and associated Aaw-Sw relationships (c-d) for natural sandy
soils. Filled symbols are the calculated GSA values for the media.

While additional measured Aaw-Sw relationships for natural soils would
be needed to validate the corrections made here to the Benson PTF, theseresults indicate that it is possible that a simple correction to the PTF
could render it more useful for estimating Aaw-Sw relationships by the
LTM for natural sandy soils containing fines. This discussion highlights
the importance of obtaining SWRCs that are most representative of thevadose zone soils in question when using the LTM approach to derive
Aaw-Sw relationships. It is additionally important to note that applying
the LTM approach is not limited to using VG-based models to charac­terize soil-water retention characteristics. Other models, including the
Brooks-Corey model (e.g., Oostrom et al., 2001), can also be used.
3.5. Comparisons with aqueous tracer-based methods
In addition to the XMT methods, air-water interfacial adsorbing
tracers have also been used to estimate Aaw. These have included both
gas-phase tracers (GT) and aqueous-phase tracers (AT). The GT methods
have been shown to provide Aaw values that are significantly greaterthan those provided by aqueous tracers, with Aaw approaching the areaof the media determined by BET measurement as Sw→ 0 (Peng andBrusseau, 2012). The XMT methods provide maximal Aaw values thatclosely match the GSA for uniform media with smooth surface areas. The

Applying these corrected equations, the resulting values of α for theVinton and Hayhook soils (drying curve) become 0.026 cm− 1 and 0.02
cm− 1, respectively. Likewise, the n values for these same soils become
3.53 and 2.91, respectively. As shown in Fig. 2, using these correctedparameters results in SWRCs and Aaw-Sw relationships that closely match
those derived from the measured VG parameters. This is particularly thecase for the Hayhook soil, where the fines fraction is more considerable.
6

J.A.K. Silva et al.

Journal of Contaminant Hydrology 247 (2022) 103984

fact that the Aaw-Sw relationships provided by the LTM closely match
those of the XMT methods for the uniform smooth and spherical media
categories implies that the LTM is a good predictor of the interfacial area
associated with natural pore-water drainage and imbibition (Leverett,
1941; Kibbey and Chen, 2012). The AT methods have also been shown toprovide Aaw values that were consistent with XMT results for smoothglass bead media (e.g., El Ouni et al., 2021), indicating both methods arequantifying the same interfacial areas (see, e.g., Fig. 3 – glass beadcomparison). However, for some model sands and natural soils, the AT
methods have provided Aaw values that are greater than those measured
by XMT, with maximal values falling between those measured by the GT
and XMT methods.
As demonstrated in Fig. 3, the XMT and LTM methods provide lower
Aaw values than the AT methods for the Accusand media and Vinton soil.
The two AT methods shown represent the results of both advecting
aqueous tracer (A-AT) methods and those determined by mass balance(MB) methods (Schaefer et al., 2000; Araujo et al., 2015). For the
Accusand (40/50 grade, d50 = 0.35 mm sand in all cases), the A-AT data
presented in Fig. 3 was derived from the published results of two
separate experiments using different interfacial-adsorbing tracers(including one in which the PFAS perfluorodecanoic acid (PFDA) wasused as a tracer) and performed under different flow and input con­centration conditions (Brusseau et al., 2015; Brusseau et al., 2020).
These A-AT experiments provide similar Aaw-Sw relationships (Fig. 3).
However, the extent to which these A-AT experiments can characterize
Aaw-Sw relationships is limited to the higher range of Sw (e.g., Sw > 0.5),
given the difficulty in maintaining a steady and uniform flow condition
at lower Sw for these types of experiments.
The MB method results derived from the published results of
Schaefer et al. (2000) generally exhibit the same linear trending Aaw-Sw
relationship for the Accusand media as that provided by the A-AT re­sults. The MB method is considered a viable approach to quantifying Aaw
at lower Sw conditions because it does not rely on flow, but rather on theequilibrium distribution of tracer between the interface and otherenvironmental compartments. However, the approach does result in
more scatter in the Aaw-Sw data. Brusseau (2020) noted that the scatter inthe MB method is a potential limitation of the method. The results
provided by Araujo et al. (2015), also presented in Fig. 3 for the same
Accusand sand, exhibit considerably more scatter in the data, with Aaw
values that are considerably greater than those provided by Schaefer
et al. (2000). While it might be tempting to plot a line through the cloudof datapoints, this approach is not tenable given that this ignores theresults of Schaefer et al. (2000), who used a similar approach and thesame sand. Results similar to Schaefer et al. (2000) were reported byChen and Kibbey (2006) for F110 Ottawa sand undergoing multiple
draining cycles using sodium octylbenzene sulfonate as a surface-active
tracer and an alternative measurement method. The uncertainty in the
MB-measured Aaw values for the Vinton sand is even more significant. As
a result, going forward, we will rely on the A-AT dataset as a basis forcomparing results.
The disparity between the XMT data and the A-AT Aaw values has
been attributed to the ability of the tracer to access the interfacial area
associated with micro-scale rough particle surfaces that is too small to beresolved by the XMT methods (e.g., Brusseau et al., 2006; Brusseau et al.,2007; McDonald et al., 2016). That is, the A-AT method is characterizing
both Aaw associated with capillary-held pendular water and some
additional fraction of the Aaw associated with pore surface films of waterthat is impacted by the microscale roughness of the surface. However,
the effect of surface roughness on the magnitude of the total Aaw existing
within an unsaturated porous media is a topic of continuing researchand debate.
For example, using stereoscopic scanning electron microscopy (SEM)
methods and surface configuration modeling, Kibbey (2013) demon­
strated that surface-associated water remains held to the surface via
capillary tension, which results in a significantly smoother air-waterinterfacial surface than that of the underlying solid surface. As aresult, the impact of surface roughness on the Aaw was shown to begreatly minimized during drainage for capillary pressure heads less than
− 100 cm. A capillary pressure head of − 100 cm or less would be neededto drain most of the media included in this evaluation to residual satu­
ration, including the Vinton soil. This would suggest that these same
media should present a smooth air-water interface for these surfaces
under normal drainage, with the magnitude of Aaw associated with thesewetted surfaces being well-characterized by the XMT methods. Further,
if water associated with surface roughness is held by capillarity, thepresence of this water should be represented in the SWRC, as demon­
strated in the literature (e.g., Dullien et al., 1989; Zheng et al., 2015).
Therefore, the LTM Aaw-Sw relationship predictions should be addi­
tionally representative. This was also the position of Kibbey and Chen
(2012), who contend that the LTM provides the maximum possible area
formed during natural drainage.
The utility of the A-AT method to characterize interfacial areas wasFig. 3. Comparison of XMT, LTM, and AT derived Aaw-Sw relationships forporous media with different degrees of surface roughness.
7

J.A.K. Silva et al.

Journal of Contaminant Hydrology 247 (2022) 103984

further challenged by Kibbey and Chen (2012), who used pore networkmodeling to demonstrate that neither diffusion nor head-driven flow
was sufficient to allow access to surface films (i.e., those presumed to be
associated with the increased Aaw for A-AT). These authors also sug­gested that, in the absence of additional experimental artifacts, surface
charge or surface tension gradient-driven advection would be needed forsurface-active tracers to access these surface films. This possibility was
examined by Costanza-Robinson and Henry (2017), who performed anumerical modeling study to demonstrate the impact of surface-tension
induced flow potentially providing artificially enhanced Aaw values
determined by the A-AT methods. However, given the reported repro­
ducibility of the A-AT results for experiments using different tracers at
different concentrations (Brusseau et al., 2015; Brusseau et al., 2020),
and the recent method validity evaluation provided by El Ouni et al.
(2021), who demonstrated that bulk surface-tension induced flow andother potential experimental artifacts were not relevant, it is difficult toignore the A-AT measurements. In their modeling evaluation, Kibbey
and Chen (2012) did not include the effects of microscopic surface
roughness, an omission that was thereafter addressed by Jiang et al.
(2019), who showed that including surface roughness resulted in thesame monotonic increases in Aaw values at high to moderate Sw values
and exponentially increasing Aaw values at low Sw, as observed by the GT
methods.
However, while Jiang et al. (2019) demonstrate that these higher Aaw
values could be caused by surface roughness at lower Sw and that theresulting Aaw-Sw relationships can be made to match measured data,
they did not extend this evaluation to include a reason as to how thesesurface-active tracers would contact this additional surface roughnessassociated interfacial area during aqueous transport under an assumed
uniform and constant soil moisture condition. Kibbey (2013) demon­
strated a thicker layer of water associated with rough surfaces andsuggested that diffusion of solutes within this thicker capillary-held
surface-associated water may be possible, wherein diffusive rates may
increase by roughly two orders of magnitude due to the resulting larger
cross-sectional areas for diffusion. This diffusive mechanism could, inturn, allow surface-active tracers to contact a larger portion of the totalinterfacial area, thus yielding larger AT-measured Aaw values thanwould be the case for a smooth surfaced media. Diffusion of solutes intoand out of these thickened surface-associated water films could havecontributed to the rate-limited diffusional mass transfer characterized in
the column experiments presented by Brusseau (2020) for the waterunsaturated transport of non-reactive and surface-active solutes be­
tween advective and non-advective domains. Therefore, it is possiblethat this thickened-film diffusive mass transfer could be contributing to
the disparity between the Aaw values provided by the A-AT and LTM.
Alternatively, and despite the convincing discussions provided by El
Ouni et al. (2021) and Brusseau et al. (2020), it seems likely that asurface-active tracer adsorbing to an air-water interface would impact
the magnitude of the interfacial area it is accesing for measurement. The
A-AT method is physically a chromatographic technique that assumes
that the stationary phase (i.e., the air-water interface in this case) has aconstant area spatially and temporally during the test. While this can bereliably assumed to be the case for solid adsorbent media, fluid-fluid
interfacial areas are generally not as stable. In the context of the pre­vious discussion, perturbations in the surface charge or surface pres­sure/tension of the air-water interface due to interfacial adsorption,even at a very low concentration, could still disrupt cohesive forces
between water molecules at the interface enough to modify the config­
uration of water films associated with rough surfaces. These perturba­
tions might not be significant enough to drain interparticle pendular
water but could cause a localized displacement of water films, ordeformation of water film surfaces, that could manifest as a localized
increase in Aaw as the pulse of the surface-active tracer moves throughthe porous medium. For the A-AT method, this increased interfacial area
would increase tracer retention and ultimately provide a larger
measured Aaw value than would be predicted by the LTM at the same Sw.

It is unlikely that this type of film water displacement would beobservable at the column or continuum scale. In this context, the A-AT
method would then provide an effective Aaw for these amphiphilic
surface-active tracers and not necessarily the true Aaw for the unsatu­
rated medium, which has been reported to be provided by the LTM(Kibbey and Chen, 2012). As shown in Fig. 3, when the surface rough­ness is absent, both the A-AT and LTM methods predict the same Aaw-Sw
relationship, which strengthens the argument that surface roughness is a
key factor in the A-AT results.
3.6. Modifying the LTM to include tracer-accessible area
While the causal mechanisms are not fully understood, the experi­mental results presented in the literature suggest that aqueous interfa­
cial tracers are measuring higher Aaw values than those predicted by the
LTM, and the discrepancy appears to be related to changes in the degreeof surface roughness of the porous media (see Fig. 3, Table 1). This is notnecessarily problematic for PFAS because these chemicals, like thetracers used for the A-AT method, are also surface-active and amphi­
philic. Therefore, Aaw-Sw relationships that are representative of surfaceactive solutes are needed for predicting the retention contribution of airwater interfacial adsorption during transport. What is problematic isthat continuum models, like the LTM, do not account for the observed
effects of surface-active solutes and surface roughness. The following
modification to the LTM is proposed to correct this discrepancy.
The modification applies a direct surface roughness multiplier (SRM)
to the Aaw-Sw relationship provided by the LTM (recall Eq. (1)). This
simple approach allows the LTM prediction to continue to include theeffects of particle size, uniformity, and differences in packing density as
it relates to particle angularity – as governed by the SWRC for capillarywater – while also capturing the contributive effects of surface rough­ness on the magnitude of Aaw at a given moisture condition. Using the AAT data as the benchmark Aaw-Sw relationship in this application, theSRM values needed to match the A-AT Aaw-Sw relationship for the
Accusand and Vinton soil examples are 2.4 and 3.3, respectively.
The results of applying the SRM to the Accusand and Vinton soil AawSw relationships are also presented in Fig. 3. Note that the SRM-adjusted
LTM result for the Accusand media maintains a more linear trending
Aaw-Sw relationship with decreasing Sw, owing to the particle size uni­
formity of the model sand. This result is consistent with the MB tracer
results of Schaefer et al. (2000), the sphericity and uniformity of thesand, and the low soil surface roughness factor (SRF) calculated for thissand (i.e., SRF = 18 for the medium Accusand as in Table 1). In contrast,
the SRM-adjusted LTM result for the natural Vinton soil exhibits a more
distinct non-linearly increasing Aaw-Sw relationship with decreasing Sw,
consistent with the increased fines fraction, less uniform particle-size
distribution, and factor of 18.6 increase in the SRF for this soil ascompared to the model sand. It is important also to note that thisapproach does not attempt to characterize the significantly increased
Aaw values observed by the GT results at lower Sw. In keeping with thelimitations of the LTM (Leverett, 1941; Bradford and Leij, 1997), the low
Sw endpoint of the SRM-adjusted Aaw-Sw relationships are truncated at
the residual saturation for a given porous media, which is appropriate
for unsaturated transport under aqueous advection (i.e., advective flow
is not expected when water content is below residual saturation).
To guide the selection of SRM values for other natural soils, a soil
surface roughness factor (SRF), initially proposed by Wenzel (1936) andapplied more recently by others when evaluating Aaw-Sw relationships(e.g., Brusseau et al., 2007; Costanza-Robinson et al., 2008; Jiang et al.,2019), was used. As additionally described in Table 1, the SRF is rede­fined here as:
SRF = N2 BET Area/GSA(11)

where GSA is the smooth geometric surface of the soil calculated using
Eq. (4). Properties for 21 Arizona soils were collected from the literature
8

J.A.K. Silva et al.

Journal of Contaminant Hydrology 247 (2022) 103984(i.e., Arthur et al., 2012; Ghanbarian et al., 2021). This data is provided
in Table 2. In addition to the SRF, Ghanbarian et al. (2021) provided
measured soil surface fractal dimensions (Ds) for each of these same soils
from water retention curves that were used to additionally characterizesoil surface roughness. Ds is an index of the topographical complexity ofsoil surfaces as it relates to scale and characterizes the change in surface
roughness between soil textural classifications. As these Ds values were
derived from water retention data, they are representative of the degreeof surface roughness at the continuum scale, in keeping with the LTM
predictions. The three-dimensional Ds values vary between 2 and 3(Pfeifer et al., 1983; Wang et al., 2005). A value of 2 indicates a smoothsurface, and a value of 3 corresponds to an extremely rough soil particle
surface. As observed in Table 2, Ds values generally increase withincreasing clay content. Ds values ranging between 2.09 and 2.3 havebeen documented for model sands similar to those described in Table 1(Avnir et al., 1985; Yang et al., 2016).
The Ds values were used to approximate the rate of change of surface
roughness for these soils relative to their corresponding SRF values (asshown in Fig. 4a). This slope was then used, in concert with the SRFs and
SRMs determined for the Accusand and Vinton soil, to provide data withwhich to develop an empirical model to estimate SRM values for soils ofvaried texture. A Langmuir-like equation was used to fit the Accusand
and Vinton data and the slope of the Ds-SRF relationship was used toscale the extrapolation of the model with increasing SRF. The results ofthis work are presented in Fig. 4b. Note that the proposed SRM-SRF
model is scaled such that SRM = 1 when SRF = 1, which defines asmooth surface condition consistent with the Aaw-Sw results for the glass
bead case shown in Fig. 3. Additionally, in keeping with the trend pro­vided by the Ds-SRF relationship, SRM values do not exceed a value of
3.5 as SRF values increase.
Unfortunately, no additional Aaw-Sw data for natural soils was foundin the published literature, precluding further validation of the proposed
SRF-SRM model. This lack of Aaw-Sw relationship data presents a sig­
nificant data gap for this area of research and is problematic with respect
to assessing the fate and transport of surface-active contaminants, like
PFAS, as there is a growing need to model PFAS transport in vadosezones comprised of soils of varied textural classification. Strictly, theSRF-SRM model would be most applicable to sands and predominantly
sandy soils similar to those described herein. On the other hand, for soils
with broad particle size distributions (e.g., loams and clayey soils), it ispossible that the effects of surface roughness on the Aaw-Sw relationshipcould be overwhelmed by the contributions of interparticle capillaryeffects, obviating the need for correcting the LTM prediction. Therefore,
more research is needed to determine the importance of surface roug­ness on Aaw values for a wider variety of natural soils and sediments,
particularly with respect to supporting site-scale modeling. However,the proposed SRF-SRM model provides a practical first-approximation
method that incorporates the effects of surface roughness on Aaw esti­mates provided by the LTM, where the SRM can be estimated from
common soil textural classification measurements (i.e., the BET specificsurface area of the soil sample and the median particle size, or d50).
Examples of LTM-predicted Aaw-Sw relationships are compared inFig. 5 to their corresponding SRM-corrected relationships for a few
sandy soils and a loam selected from Table 2. Textural data from Table 2
was used as input to the Rosetta model (described previously) to provide
VG parameters and generate SWRCs that were then used to generate
LTM predictions. Texturally, the sand and loamy sand results are mostconsistent with the Vinton soil and the Hayhook soil, as are the predicted
Aaw-Sw relationships. Regarding the loamy sand predictions, the SRMadjusted Aaw-Sw relationship closely matches that of the SRM-adjusted
sand results until Sw= 0.5. Aaw values increase considerably thereafter,
representing the effect of surface roughness on the magnitude of the Aaw
with increased pore-water drainage at the continuum scale. As the fines
content of the soil increase, predicted Aaw values increase. For the sandy
loam and loam soils, predicted Aawvalues are roughly a factor of 10
greater than those of the sand and loamy sand as Sw falls below 0.2. Inthese cases, it is important to consider that field capacities for a typicalsandy loam and loam soil range between θ = 0.18–0.28 and θ =0.25–0.32, respectively (or Sw = 0.47–0.73 and Sw = 0.6–0.76, respec­tively). Under typical field conditions, these soils should retain enoughmoisture to limit the magnitude of the Aaw and the contribution of AWI
adsorption during PFAS transport. However, this would not necessarily
hold true for soils positioned at the ground surface in an arid environ­ment, where much lower water saturations could be achieved via
evaporative losses, where Aaw values could be large, and PFAS could besignificantly retained within the shallow vadose zone as a result.
3.7. The significance of the Aaw-Sw relationship on transport
The potential effect of this Aawvariability on PFAS retardation factors(R) is demonstrated for PFOS, a representative PFAS, and is presented inFig. 6. R-values were calculated using the following equation:

Table 2
Soil properties data used to modify LTM prediction.
Arizona soils

SSABETa (cm− 1)

Category%sand/silt/clay

porosityd50 (mm)

GSA (cm− 1)

SRF

SRM

Ds (− )

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

22,80016,80037,50093,75076,500111,000
296,100
107,550
149,400
150,750
222,900
466,800
392,550
317,700
525,600
436,050
397,050535,200
781,500
570,750
92,250sand
sand
loamy sandsandy loam
loamy sandsandy loam
sandy loam
sandy loam
sandy clay loam
silt loam
loam
loam
sandy clay loam
sandy clay loam
sandy clay loam
silt loam
clay loam
silty clay
clay
sandy clay loam
sandy loam

91/7/2
96/3/181/14/5
74/18/8
80/14/6
67/25/8
71/16/13
59/32/9
53/26/21
22/53/25
37/46/17
39/40/21
66/15/27
58/15/27
51/22/27
4/73/23
26/45/29
9/40/51
29/19/52
49/21/30
71/20/9

0.36
0.36
0.38
0.38
0.37
0.39
0.39
0.4
0.43
0.46
0.42
0.4
0.4
0.42
0.43
0.45
0.43
0.43
0.45
0.41
0.38

0.17
0.17
0.135
0.098
0.135
0.098
0.06
0.098
0.07
0.027
0.035
0.035
0.035
0.027
0.027
0.027
0.018
0.006
0.006
0.027
0.098

22.6
22.6
27.6
38.0
28.0
37.3
61.0
36.7
48.9
120.0
99.4
102.9
102.9
128.9
126.7
122.2
190.0
570.0
550.0
131.1
38.0

1009
744
1361
2470
2732
2972
4854
2928
3058
1256
2242
4538
3816
2465
4149
3568
2090
939
1421
4353
2430

3.43
3.41
3.45
3.47
3.47
3.48
4.31
4.11
4.13
3.78
4.01
4.28
4.21
4.04
4.25
4.19
3.98
3.67
3.83
4.27
4.04

2.65
2.68
2.71
2.68
2.70
2.71
2.76
2.69
2.74
2.74
2.70
2.76
2.78
2.81
2.81
2.76
2.78
2.81
2.87
2.83
2.70

a

Reported specific surface area from N2/BET. d50 values from Dosskey et al. (2006). Ds is the surface fractal dimension(Ghanbarian et al., 2021).
9

J.A.K. Silva et al.

Journal of Contaminant Hydrology 247 (2022) 103984

Fig. 4. Data used to develop an empirical model to estimate SRM values: (a) the rate of change of Ds with increasing SRF, and (b) model fit to available SRM-SRF
data. Dashed line in 4b indicates extrapolation of model prediction to higher SRF based on the slope of the surface fractal dimension in 4a.

Fig. 5. LTM and SRM-adjusted LTM predictions of the Aaw-Sw relationships for media selected from Table 2: (a) sand 2 and loamy sand 5, (b) sandy loam 6 andloam 11.

Fig. 6. Calculated ranges for PFOS retardation with changing Sw for the Accusand media and the Vinton soil. Combined refers to the inclusion of solid-phase sorption
and AWI adsorption. LTM refers to the use of the LTM model (Eq. (1)) and SRM is the LTM model result when the SRM is applied.

10

J.A.K. Silva et al.

ρ Kd kaw Aaw
R = 1+ b +θθJournal of Contaminant Hydrology 247 (2022) 103984

media become rougher (as characterized by the proposed SRF). The AAT method was shown to provide Aaw values that are consistently
greater than both the XMT and LTM results for a given Sw, with thedegree of difference increasing with increased surface roughness,consistent with the conclusions drawn from literature sources whereinmuch of the A-AT data was obtained.
The A-AT method utilizes surface-active solutes as tracers to char­acterize Aaw at the continuum scale in porous media columns. These
tracers appear to be sampling portions of the air-water interface that are
not represented by either the LTM prediction or XLM measurements.
While the exact mechanism for this disparity is uncertain, it appears thatthe degree of microscale surface roughness is a dominant contributing
factor. Because the AT method uses tracers that are surface-active, andthe contaminants of interest in this work are also surface-active, it isreasonable to conclude the tracers themselves may influence the inter­
facial area. Thus, the A-AT data provides the most representative Aawvalues for surface-active solutes like PFAS. In addition, the LTM (andsimilar thermodynamic methods) likely provide a true measure of the
Aaw under natural drainage conditions, a conclusion based on the pre­sentations of Kibbey and Chen (2012) and Kibbey (2013) and on thesimilarity in predicted Aaw-Sw relationships when media surfaces aresmooth.
Using the literature-derived A-AT measured Aaw-Sw relationships as abenchmark for surface-active tracers, the LTM prediction was modified
to include the effects of surface roughness on the predicted Aaw-Sw
relationship by applying a direct multiplier to the LTM prediction toaccount for surface roughness (i.e., the SRM). An empirical model wasdeveloped using the measured data to estimate the magnitude of the
SRM with changes in calculated SRF values as soil textural classifications
change. Because this empirical model was developed from a limited
dataset, it should be considered to provide only first-order estimates;
additional measured Aaw-Sw relationships for a wider variety of soil
textures are needed for a more robust test of this empirical model.
However, in the absence of this additional data, the approach provides apractical method for calculating or simulating PFAS mass transport in
the vadose zone and accounting for the effects of surface roughness on
Aaw-Sw relationships.
Finally, because the representativeness of the Aaw-Sw relationships
predicted by the LTM depends on the representativeness of the SWRC
used, methods of estimating VG parameters that define the SWRCs werecompared. The accuracy of the Rosetta model was found to improvewhen fines (e.g., particles <0.0075 mm) are present in the soil. For nonuniform sandy soils without fines, the Benson PTF, which corrects pre­dictions for changes in the uniformity of porous media particles, wasshown to best match measured SWRCs for the sands presented in thiswork. The Benson PTF was also modified to account for the presence offines, which greatly improved the accuracy of the Aaw-Sw relationshipprediction for the natural soils used in this work compared to thatderived from the measured SWRCs. However, as with all PTFs, including
the Rosetta model, it is important to understand the operational limi­
tations of these models to ensure representative prediction of hydraulicproperties.

(12)

where ρb is the soil bulk density (M/L3), Kd is the solid-phase sorption
coefficient (L3/M1), kaw is the air-water adsorption coefficient (L), and θ= Swn. As presented, Eq. (12) assumes linear solid-phase sorption andthat sorption and air-water interfacial adsorption are the only sources ofretention. Sorption coefficients (Kd) were 0.16 and 1.2 cm3/g for the
Accusand and Vinton soil, respectively, and were derived from the linear
model fit to the data presented by van Glubt et al. (2021) for soluteconcentrations <1 mg/L. Other parameters used in these calculations
were: ρb= 1.5 g/cm3 and kaw= 0.026 cm. The results of these calcula­tions (Fig. 6) demonstrate not only the variability in PFOS retention, butalso the importance of using the appropriate Aaw-Sw relationship whensurface-active solutes are being considered, because vadose-zone PFAS
retention could be greatly misrepresented.
Finally, it is important to note that all Aaw-Sw relationships used inthis work were derived from data representing drainage conditions. The
magnitude of Aaw values in the vadose zone will be impacted by tran­sient cycles of drainage and imbibition (wetting) during natural infil­
tration events. While Leverett (1941) demonstrates the LTM could beapplied to both drainage and imbibition conditions, characterized byseparate SWRCs and potentially providing different Aaw-Sw relationships
for each condition, the validity of this application with respect to thetransport of PFAS has not been. Using A-AT methods and the Vinton soil,
Brusseau et al. (2007) observed that Aaw-Sw relationships measured
under both drainage and imbibition conditions were similar and citedthe importance of surface roughness associated water films for thesimilarity in measured Aaw values. Schaefer et al. (2000) also observed
similar Aaw-Sw relationships under drainage and imbibition conditions
using MB methods. However, both studies utilized a surface-active
tracer (i.e., sodium dodecyl benzenesulfonate) at concentrations thatwould promote interfacial tension reduction and could have altered the
moisture retention characteristics that would have been observed forwater without surfactant. Therefore, additional research is needed on
this topic.
4. Conclusions
When performing calculations or numerical simulations of the fate
and transport of PFAS and other surface-active solutes in the vadosezone, and accurate Aaw-Sw relationship is needed. Most research has
focused on accurate characterization of interfacial adsorption co­efficients and the associated concentration dependence. However,
capturing changes in the Aaw-Sw relationship relationship that result
from changes in soil texture is equally important. Indeed, accurate
determination of the Aaw-Sw relationship is arguably more importantwhen one considers the challenges in characterizing heterogeneity atboth the pore- and continuum scale and its relationship to the SWRC and
Aaw-Sw relationship. In this work, the well-known LTM model iscompared to two methods for the determination of the Aaw-Sw rela­tionship: microtomography and aqueous tracer methods. The LTM
model is favored here due to its simplicity, and because it represents thechange in the magnitude of the Aaw with changes in Sw at the continuum
scale, consistent with practical methods used to model PFAS transport (i.
e., continuum mechanics behind the SWRC and the application of the
Richards equation for unsaturated flow). In addition, the LTM method isless time-consuming and expensive, yet appropriately captures the ef­
fects of soil particle size, uniformity in the particle size distribution, anddifferences in packing density on the Aaw-Sw relationship.
The XMT, the A-AT method, and the LTM were found to provide very
similar Aaw-Sw relationships when the porous media surface is smooth (i.
e., glass beads), suggesting each method characterizes the same effective
interfacial area in these cases. However, both the LTM and A-AT
methods provide larger Aaw values than measured by the XMT method
for non-uniform, increasingly angular media, or when the surfaces of theCRediT authorship contribution statement
Jeff A.K. Silva: Conceptualization, Methodology, Validation, Formal
analysis, Investigation, Resources, Writing – original draft, Writing –review & editing, Visualization, Supervision, Project administration,
Funding acquisition. Jiří Šimůnek: Methodology, Software, Resources,
Writing – review & editing. John E. McCray: Methodology, Writing –review & editing.
Declaration of Competing Interest
The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence11

J.A.K. Silva et al.

Journal of Contaminant Hydrology 247 (2022) 103984

the work reported in this paper.

Ghanbarian, B., Hunt, A.G., Bittelli, M., Tuller, M., Arthur, E., 2021. Estimating specificsurface area: Incorporating the effect of surface roughness and probing moleculesize. Soil. Sci. Am. J. https://doi.org/10.1002/saj2.20231. Advance online
publication.
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. https://doi.org/10.1029/2007WR006038.
Guelfo, J.L., Higgins, C.P., 2013. Subsurface transport potential of perfluoroalkyl acids ataqueous film-forming foam (AFFF)-impacted sites. Environ. Sci. Technol. 47,
4164–4171.
Guo, B., Zeng, J., Brusseau, M.L., 2020. A mathematical model for the release, transport,
and retention of per- and polyfluoroalkyl substances (PFAS) in the vadose zone.
Water Resour. Res. 56 (2), e2019WR026667 https://doi.org/10.1029/
2019WR026667.
Gvirtzman, H., Roberts, P.V., 1991. Pore scale spatial analysis of two immiscible fluids inporous media. Water Resour. Res. 27 (6), 1165–1176.
Hamamoto, S., Moldrup, P., Kawamoto, K., Sakaki, T., Nishimura, T., Komatsu, T., 2016.
Pore network structure linked by X-ray CT to particle characteristics and transportparameters. Soils Found. 56 (4), 676–690.
Jiang, H., Guo, B., Brusseau, M.L., 2019. Pore-scale modeling of fluid-fluid interfacial
area in variably saturated porous media containing microscale surface roughness.
Water Resour. Res. 56, e2019WR025876 https://doi.org/10.1029/2019WR025876.
Kakare, M.V., Fort, T., 1996. Determination of the air-water interfacial area in wet“unsaturated” porous media. Langmuir 12 (8), 2041–2044.
Kibbey, T.C.G., 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. https://doi.org/10.1002/wrcr.20383.
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. https://doi.
org/10.1029/2012WR011862.
Kim, H., Rao, P.S.C., Annable, M.D., 1997. Determination of effective air-waterinterfacial area in partially saturated porous media using surfactant adsorption.
Water Resour. Res. 33 (12), 2705–2711.
Kim, H., Rao, P.S.C., Annable, M.D., 1999. Gaseous tracer technique for estimating airwater interfacial areas and interface mobility. Soil Soc. Am. J. 63, 1554–1560.
Leverett, M.C., 1941. Capillary behavior in porous solids. Trans. Am. Inst. Min. Merall.
Pet. Eng. 142, 152–169.
McDonald, K., Carroll, K.C., Brusseau, M.L., 2016. Comparison of fluid-fluid interfacial
areas measured with x-ray microtomography and interfacial partitioning tracer tests
for the same samples. Water Resour. Res. 52 (7), 5393–5399. https://doi.org/
10.1002/2016WR018775.
McKenzie, E.R., Siegrist, R.L., McCray, J.E., Higgins, C.P., 2016. The influence of a
nonaqueous phase liquid (NAPL) and chemical oxidant application of perfluoroalkyl
acid (PFAA) fate and transport. Water Res. 92, 199–207.
Oostrom, M., White, M.D., Brusseau, M.L., 2001. Theoretical estimation of free andentrapped non-wetting fluid interfacial areas in porous media. Adv. Water Resour.
24, 887–898.
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.
Peng, S., Brusseau, M.L., 2005. Impact of soil texture on air-water interfacial areas inunsaturated sandy porous media. Water Resour. Res. 41, W03021. https://doi.org/
10.1029/2004WR003233.
Peng, S., Brusseau, M.L., 2012. Air-water interfacial area and capillary pressure: porousmedium texture effects and an empirical function. J. Hydrol. Eng. 17 (7) https://doi.
org/10.1061/(ASCE)HE.1943-5584.0000515.
Pfeifer, P., Avnir, D., Farin, D., 1983. Ideally irregular surfaces, of dimension greater thantwo, in theory and practice. Surf. Sci. 126, 569–572. https://doi.org/10.1016/00396028(83)90759-8.
Sakaki, T., Illangasekare, T.H., 2007. Comparison of height-averaged and pointmeasured capillary pressure-saturation relations for sands using a modified Tempecell. Water Resour. Res. 43, W12502. https://doi.org/10.1029/2006WR005814.
Schaap, M.G., Leij, F.J., van Genuchten, M.Th., 2001. Rosetta: a computer program forestimating soil hydraulic parameters with hierarchical pedotransfer functions.
J. Hydrol. 251, 163–176.
Schaefer, C.E., DiCarlo, D.A., Blunt, M.J., 2000. Experimental measurement of air-waterinterfacial area during gravity drainage and secondary imbibition in porous media.
Water Resour. Res. 36 (4), 885–890.
Schroth, M.H., Ahearn, S.J., Selker, J.S., Istok, J.D., 1996. Characterization of Millersimilar silica sands for laboratory hydrologic studies. Soil Sci. Soc. Am. J. 60,
1331–1339.
Schroth, M.H., Oostrom, M., Dobson, R., Zeyer, J., 2008. Thermodynamic model for
fluid-fluid interfacial areas in porous media for arbitrary drainage-imbibition
sequences. Vadose Zone J. 7, 966–971. https://doi.org/10.2136/vzj2007.0185.
Shin, H.-M., Viery, V.M., Ryan, P.B., Detwiler, R., Sanders, B., Steenland, K., Bartell, S.
M., 2011. Environmental fate and transport modeling for perfluorooctanoic acidemitted from the Washington works Facility in West Virginia. Environ. Sci. Technol.
45 (4), 1435–1442.
Silva, J.A.K., 1997. Retention Processes Affecting VOC Vapor Transport in WaterUnsaturated Porous Media. MSc Thesis. University of Arizona, Tucson, AZ.
Silva, J.A.K., Bruant, R.G., Conklin, M.H., Corley, T.L., 2002. Equilibrium partitioning ofchlorinated solvents in the vadose zone: low foc geomedia. Environ. Sci. Technol. 36(7), 1613–1619.

Acknowledgements
This material is based on work supoored by the U.S. Army Corps ofEngineers and The Department of Defense Strategic Environmental
Research and Development Program (SERDP) under Contract No.
W912HQ18C0076.
References
Anderson, R.H., Long, G.C., Porter, R.C., Anderson, J.K., 2006. Occurrence of selectperfluoroalkyl substances at U.S. Air force aqueous film-forming foam release sites
other than fire252 training areas: field-validation of critical fate and transport
properties. Chemosphere 150, 678–685. https://doi.org/10.1016/j.
chemosphere.2016.01.014.
Anderson, R.H., Adamson, D.T., Stroo, H.F., 2019. Partitioning of poly- andperfluoroalkyl substances from soil to groundwater within aqueous film-forming
foam source zones. J. Contam. Hydrol. 220, 59–65. https://doi.org/10.1016/j.
jconhyd.2018.11.011.
Araujo, J.B., Brusseau, M.L., 2020. Assessing XMT-measurement variability of air-waterinterfacial areas in natural porous media. Water Resour. Res. 56, e2019WR025470
https://doi.org/10.1029/2019WR025470.
Araujo, J.B., Mainhagu, J., Brusseau, M.L., 2015. Measuring air–water interfacial area forsoils using the mass balance surfactant-tracer method. Chemosphere 134, 199–202.
Arthur, E., Tuller, M., Moldrup, P., Resurreccion, A.C., Meding, M.S., Kawamoto, K.,
Komatsu, T., de Jonge, L.W., 2012. Soil specific surface area and non-singularity ofsoil-water retention at low saturations. Soil. Sci. Am. J. 77, 43–53.
Avnir, D., Farin, D., Pfeifer, P., 1985. Surface geometric irregularity of particulatematerials: the fractal approach. J. Colloid Interface Sci. 103 (1), 112–123.
Benson, C.H., Chiang, I., Chalermyanont, T., Sawangsuriya, A., 2014. Estimating van
Genuchten parameters α and n for clean sands from particle size distribution data. In:
Geotechnical Special Publication No. 233, From Soil Behavior Fundamental toInnovations in Geotechnical Engineering. American Society of Civil Engineers.
Bradford, S.A., Leij, F.J., 1997. Estimating interfacial areas for multi-fluid soil systems.
J. Contam. Hydrol. 27, 83–106.
Bradford, S.A., Wang, Y., Torkzaban, S., Šimůnek, J., 2015. Modeling the release of E.
coli D21g with transients in water content. Water Resour. Res. 51, 3303–3316.
Brusseau, M.L., 2018. Assessing the potential contributions of additional retention
processes to PFAS retardation in the subsurface. Sci. Total Environ. 613-614,
176–185.
Brusseau, M.L., 2020. Simulating PFAS transport influenced by rate-limited multiprocess retention. Water Res. 168, 115179.
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.
Brusseau, M.L., Peng, S., Schnaar, G., Murao, A., 2007. Measuring air-water interfacial
areas with x-ray microtomography and interfacial partitioning tracer tests. Environ.
Sci. Technol. 41 (6), 1956–1961.
Brusseau, M.L., El Ouni, A., Araujo, J.B., Zhong, H., 2015. Novel methods for measuring
air–water interfacial area in unsaturated porous media. Chemosphere 127, 208–213.
Brusseau, M.L., Yan, N., Van Glubt, S., Wang, Y., Chen, W., Lyu, Y., et al., 2019.
Comprehensive retention model for PFAS transport in subsurface systems. Water
Res. 148, 41–50.
Brusseau, M.L., Lyu, Y., Guo, B., 2020. Low-concentration tracer tests to measure airwater interfacial area in porous media. Chemosphere 250, 126305.
Chen, L., Kibbey, T.C.G., 2006. Measurement of air-water interfacial area for multiple
hysteretic drainage curves in an unsaturated fine sand. Langmuir 22 (16),6874–6880. https://doi.org/10.1021/la053521e.
Costanza-Robinson, M.S., 2001. Elucidation of Retention Processes Governing theTransport of Volatile Organic Compounds in Unsaturated Soil Systems. Ph.D.
Dissertation. University of Arizona.
Costanza-Robinson, M.S., Henry, E.J., 2017. Surfactant-induced flow compromises
determination of air-water interfacial areas by surfactant miscible-displacement.
Chemosphere 171, 275–283.
Costanza-Robinson, M.S., Harrold, K.H., Lieb-Lappen, R.M., 2008. X-ray
microtomography determination of air-water interfacial area-water saturation
relationships in sandy porous media. Environ. Sci. Technol. 42 (8), 2949–2956.
Costanza-Robinson, M.S., Estabrook, B.D., Fouhey, D.F., 2011. Representative
elementary volume estimation for porosity, moisture saturation, and air-waterinterfacial areas in unsaturated porous media: data quality implications. Water
Resour. Res. 47, W07513. https://doi.org/10.1029/2010WR009655.
Dalla, E., Hilpert, M., Miller, C.T., 2002. Computation of the interfacial area for two-fluid
porous medium systems. J. Contam. Hydrol. 56 (1–2), 25–48.
Dosskey, M.G., Helmers, M.J., Eisenhauer, D.E., 2006. An approach for using soil surveys
to guide the placement of water quality buffers. J. Soil Water Conserv. 61 (6),334–354.
Dullien, F.A., Zarcone, C., MacDonald, I.F., Collins, A., Bochard, R.D.E., 1989. The effectsof surface roughness on the capillary pressure curves and the heights of capillary risein glass bead packs. J. Colloid Interface Sci. 127 (2), 362–372.
El Ouni, A., Guo, B., Zhone, H., Brusseau, M.L., 2021. Testing the validity of the miscibledisplacement interfacial tracer method for measuring air-water interfacial area:
independent benchmarking and mathematical modeling. Chemosphere 263, 128193.

12

J.A.K. Silva et al.

Journal of Contaminant Hydrology 247 (2022) 103984
Wang, K., Zhang, R., Wang, F., 2005. Testing of pore-solid fractal model for the soil water
retention function. Soil Sci. Soc. Am. J. 69, 776. https://doi.org/10.2136/
sssaj2004.0247.
Wang, J.-P., Lambert, P., De Kock, T., Cnudde, V., Francois, B., 2019. Investigation of theeffect of specific interfacial area on strength of unsaturated granular materials by xray tomography. Acta Geotech. 14, 1545–1559.
Weber, A.K., Barber, L.B., LeBlanc, D.R., Sunderland, E.M., Vecitis, C.D., 2017.
Geochemical and hydrologic factors controlling subsurface transport of poly- andperfluoroalkyl substances, Cape Cod, Massachusetts. Environ. Sci. Technol. 51,
4269–4279.
Wenzel, R.N., 1936. Resistance of solid surfaces to wetting by water. J. Ind. Eng. Chem.
28 (8), 988–994.
Xiao, F., Simcik, M.F., Halbach, T.R., Gulliver, J.S., 2015. Perfluorooctane sulfonate(PFOS) and perfluorooctanoate (PFOA) in soils and groundwater of a US
metropolitan area: migration and implications for human exposure. Water Res. 72,
64–74.
Yang, H., Rahardjo, H., Leong, E., Fredlund, D., 1994. Factors affecting drying andwetting soil-water characteristic curves of sandy soils. Can. Geotech. J. 41, 908–920.
Yang, H., Baudet, B.A., Ya, T., 2016. Characterization of the surface roughness of sandparticles using an advanced fractal approach. Proc. R. Soc. A 472, 20160524.
Zheng, W., Yu, X., Jin, Y., 2015. Considering surface roughness effects in a triangularpore space model for unsaturated hydraulic conductivity. Vadose Zone J. 14 (7)https://doi.org/10.2136/vzj2014.09.0121 vzj2014.09.0121.

Silva, J.A.K., Martin, W.A., Johnson, J., McCray, J.E., 2019. Evaluating air-water and
NAPL-water interfacial adsorption and retention of Perfluorocarboxylic acids within
the vadose zone. J. Contam. Hydrol. 223, 103472.
Silva, J.A.K., Šimůnek, J., McCray, J.E., 2020. A modified HYDRUS model for simulating
PFAS transport in the vadose zone. Water 12, 2758. https://doi.org/10.3390/
w12102758.
Šimůnek, J., van Genuchten, M.Th., Šejna, M., 2016. Recent developments andapplications of the HYDRUS computer software packages. Vadose Zone J. 15 (7), 25.
https://doi.org/10.2136/vzj2016.04.0033.
Sweijen, T., Aslannejad, H., Hassanizadeh, M., 2017. Capillary pressure-saturation
relationships for porous granular materials: pore morphology method vs. pore unit
assembly method. Adv. Water Resour. 107, 22–31.
Torrens, J.L., Herman, D.C., Miller-Maier, R., 1998. Biosurfactant (Rhamnolipid)
sorption and the impact on Rhamnolipid-facilitated removal of cadmium fromvarious soils under saturated flow conditions. Environ. Sci. Technol. 32 (6),776–781.
van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.
van Glubt, S., Brusseau, M.L., Ni, Y., Huang, D., Khan, N., Carroll, K.C., 2021. Column
versus batch methods for measuring PFOS and PFOA sorption to geomedia. Environ.
Pollut. 268, 115917.

13