Journal of Environmental Management 237 (2019) 576–584

Contents lists available at ScienceDirect

Journal of Environmental Management
journal homepage: www.elsevier.com/locate/jenvman

Research article

Determination of speciﬁc LNAPL volumes in soils having a multimodal poresize distribution

T

Miguel A. Alfaro Sotoa,∗, Robert Lenhardb, Hung K. Changa, Martinus Th van Genuchtena,ca

Department of Applied Geology and Center for Environmental Studies, São Paulo State University, UNESP, Rio Claro, SP 13506-900, Brazil
13210 Palatine Hill, San Antonio, TX, USAc
Department of Nuclear Engineering, Federal University of Rio de Janeiro, UFRJ, Rio de Janeiro, RJ, Brazilb

A R T I C LE I N FO

A B S T R A C T

Keywords:
Multimodal pore systems
Free LNAPL
Recoverable LNAPLLDRMIn this paper we present modiﬁcations to previously published models for determining the speciﬁc volume ofnon-aqueous phase liquids (LNAPLs) in the subsurface at and near the groundwater table following a spill or leakfrom the soil surface. The modiﬁcations account for porous media having multimodal pore-size distributions as isoften the case with tropical soils. Data from the literature are used to show that the use of multimodal pore-size
distributions can lead to signiﬁcantly diﬀerent subsurface LNAPL speciﬁc volume predictions and possible
LNAPL recovery rates, compared to when only unimodal pore-size distributions are considered. Diﬀerences of up
to 200% are possible when the dual-porosity nature of the pore system is ignored, which can yield erroneous
estimates of the time needed to remediate LNAPLs from contaminated areas when conventional systems are
employed.

1. Introduction
Soil and groundwater systems in many areas are contaminated bylight non-aqueous phase liquids (LNAPLs) from accidental spills on thesoil surface or by leaking underground containers (e.g., Charbeneau
et al., 2000; Hernández-Espriú et al., 2012; Lenhard et al., 2018;
Sookhak Lari et al., 2018). Governmental regulatory agencies commonly require the removal, as much as feasible, of the mobile (freephase) component of a LNAPL. Consequently, the LNAPL volume in thesubsurface is an important parameter when developing potential remediation strategies and obtaining cost estimates (Wickramanayaque
et al., 1991).
Various methods for determining the subsurface volume of LNAPLs
through inspection of monitoring wells can be found in the literature.
Investigators typically use for this purpose the LNAPL thickness ﬂoating
on top of the water table inside a monitoring well. Most approaches
initially assumed that the LNAPL thickness in the subsurface was equal
to that in the monitoring well. Subsequent investigations showed that
the LNAPL thickness in a well is not a direct measurement of the LNAPL
thickness in the subsurface, but may be 2 to 10 times larger (van Dam,
1967; De Pastrovich et al., 1979; Mercer and Cohen, 1990). Several
studies (e.g., Hampton and Miller, 1988) later concluded that simple
ratios do not provide suitable predictions of the actual LNAPL volume
in the subsurface.
∗

The saturation distribution of water, LNAPL and air in the subsurface is governed by the air-LNAPL and LNAPL-water capillary pressures(Lenhard and Parker, 1990; CL.AIRE, 2014; Lenhard et al., 2017) wherewater is the wetting ﬂuid, air the nonwetting ﬂuid, and LNAPL has
intermediate wettability (i.e., water occupies the smallest pores, air thelargest pores, and LNAPL the intermediate-sized pores). Hence, subsurface LNAPL volumes depend upon the capillary pressures and thesize distribution of pores in the subsurface. Farr et al. (1990) andLenhard and Parker (1990) were ﬁrst to develop methods for evaluating
the volume of a LNAPL in a vertical slice (the LNAPL speciﬁc volume) ofwater-wet porous media using capillary pressures, ﬂuid and porous
media properties, and the LNAPL thickness in a well. Their approachassumes that ﬂuids in the subsurface are not moving vertically and that
good contact exists between the ﬂuids in the formation and the well.
These authors, among others (CL.AIRE, 2014; ITRC, 2016) showed that
the LNAPL is variably saturated above the water-saturated zone andthat LNAPL-saturated ‘pancakes’ do not exist (Newell et al., 1995).
Studies such as Wickramanayake et al. (1991) showed that the conceptual model of Farr et al. (1990) and Lenhard and Parker (1990) was
more appropriate for estimating LNAPL volumes than earlier methods(e.g., De Pastrovich et al., 1979; Hall et al., 1984).
One frequently used analytical method to predict LNAPL speciﬁcvolume is the LDRM (LNAPL Distribution and Recovery Model) software of the American Petroleum Institute (API, 2007), developed by

Corresponding author. São Paulo State University, 1515 Av. 24A, Bela Vista, Rio Claro, SP 135050-060, Brazil.
E-mail address: alfaro@rc.unesp.br (M.A. Alfaro Soto).

https://doi.org/10.1016/j.jenvman.2019.02.077
Received 26 September 2018; Received in revised form 16 January 2019; Accepted 16 February 20190301-4797/ © 2019 Elsevier Ltd. All rights reserved.

Journal of Environmental Management 237 (2019) 576–584

M.A. Alfaro Soto, et al.

formations. To conduct such a test, we ﬁrst analyzed experimental airwater data using the traditional VG equation. We then used the airwater VG parameters to estimate air-NAPL capillary pressure curves
following scaling procedures as outlined in API (2007). Results aresubsequently compared with direct air-NAPL capillary pressure measurements or curves. We note that LNAPL-water systems were nottested since we assumed that the results would be similar to those of airLNAPL systems. Finally, the scaling procedure for multimodal pore
distributions will be used in an adaptation of the API model (2007) toobtain total and recoverable volumes of LNAPLs in several soils.

Charbeneau (2003), which employs the concepts of Farr et al. (1990)
and Lenhard and Parker (1990). This software makes it possible to simulate the ability of various remediation techniques to recover freephase (mobile) LNAPLs from the subsurface. The software allows estimates of the distribution of the LNAPL in the porous medium, including
rates, volumes and times for possible LNAPL recovery, consistent withearlier studies by Parker and Lenhard (1989). However, the LDRM
model involves a series of simpliﬁcations, such as assuming homogeneity of the formation, the ﬂuids, and overall chemical conditions,while also considering equilibrium of vertical forces and ignoring hysteresis. The model hence has certain limitations in terms of predicting
the total speciﬁc LNAPL volume present and how much can be recovered.
The models by Farr et al. (1990), Lenhard and Parker (1990),
Lenhard et al. (2004), API (2007) and Jeong and Charbeneau (2014),
among others, utilize parameters of frequently used capillary pressuresaturation functions to predict LNAPL speciﬁc volumes. The two mostcommon functions are those by Brooks and Corey (1964) and van
Genuchten (1980), denoted here as the BC and VG functions, respectively. The BC and VG models apply mostly to porous media having
normal or lognormal pore-size distributions. When the BC parameters
are used in the Farr et al. (1990), Lenhard and Parker (1990), and API(2007) models, closed-form analytical equations can be derived for the
LNAPL volumes. However, when the VG model is used, numerical integration is needed to predict LNAPL volumes.
Many porous media do not have normal or lognormal pore-size
distributions, but multimodal distributions (Peters and Klavetter, 1988;
Durner, 1994; Romano et al., 2011; Wijaya and Leong, 2016; Madi
et al., 2018). Heterogeneous pore-size distributions may be a result of
the speciﬁc particle size distribution, or the creation of secondary
porosity as a result of genetic or other processes, such as physical orchemical aggregation, biological processes, and/or the eﬀects of glaciation (morain) or soliﬂuction. Such processes are not limited to unconsolidated materials, since bimodal or multimodel distributions canalso be found in structured soils and fractured rocks with large pores(fractures) interspersed with smaller, less permeable micropores of thesoil or rock matrix (Peters and Klavetter, 1988; Gerke and van
Genuchten, 1993). Many tropical soils (such as Ultisols and Oxisols)and volcanic soils (Andisols) frequently exhibit also heterogeneous
pore-size distributions (Miguel et al., 2006; Spohrer et al., 2006; Alfaro
Soto et al., 2008; Rudiyanto et al., 2013).
Mathematical models derived for normal or lognormal pore-size
distributions may not be appropriate for predicting free LNAPL distributions for multimodel pore structures. The purpose of this study is
to allow extension of the models by Farr et al. (1990), Lenhard andParker (1990), and API (2007) to porous media having multimodel
pore-size distributions, and to determine the total and recoverable volumes of LNAPL in such systems. Results obtained with the revised
models will be compared to those obtained with the classical unimodal
formulations. In this paper, we discuss the literature and data used forthis study, summarize the equations describing the capillary pressureand LNAPL equilibrium distributions in multimodal systems, and present our results.

2.1. Air-water and air-NAPL capillary pressure curves
Hanging water (liquid) column methods (ASTM, 2003a) were used
to measure capillary pressure-saturation relationship of the air-waterand air-LNAPL systems. Initial conditions for both ﬂuid systems were
full liquid saturation. Thereafter, increments of up to 50 cm capillarypressures were applied to reduce liquid saturations. For our measurements, we used aviation kerosene as the LNAPL. Air-water and airLNAPL interfacial tensions were determined using the du Nouy ring
method (model K9-Mk1, Krüss GmbH, Germany), which consists ofmeasuring the tension required to displace a platinum ring from the airwater and air-LNAPL systems under standard conditions according to
ASTM (2014).
We additionally used the ﬁlter paper method (ASTM, 2003b) tocomplement the air-water capillary pressure curves. The method consists of placing a soil sample in contact with standardized ﬁlter paper(usually Whatman No 42 or Schleincher and Shuell No 589) in a closed
atmosphere until a balance of potentials is reached. After this, themoisture content of the soil is related to its capillary pressure that is
obtained using calibration equations based on the equilibrium moisture
content of the ﬁlter paper (Campbell and Gee, 1986; ASTM, 2003b).
Soil samples used in the measurements were from soils very prevalent in the Rio Claro formation of the Brazilian state of São Paulo.
According to Graciosa et al. (2008) and Carducci et al. (2013), among
others, these soils have multimodal pore-size distributions due tostructuring processes typical of lateritic soils of Brazil.
2.2. LNAPL distributions based on inspection of monitoring wells
The relationship between the vertical distributions of water, LNAPL
and air in a porous medium and the corresponding distribution of theseﬂuids in monitoring wells can be determined provided good contact
exists between the ﬂuids in the porous medium and the wells. The relationship can be relatively complex for transient conditions (Farr et al.,1990). However, when ﬂuids are vertically static and in equilibriumwith ﬂuids in a monitoring well, then the amount of spilled LNAPL in
the formation can be estimated from the LNAPL in the well. The variables involved in this equilibrium between the ﬂuids in the well and theformation are summarized in Fig. 1a.
According to Parker and Lenhard (1989), Farr et al. (1990) andLenhard and Parker (1990), the elevations of liquids in a monitoring
well at static conditions can be used to predict capillary pressures in thesubsurface, thus making it possible to calculate the ﬂuid saturations inthe subsurface. Fig. 1b illustrates the relationship between pressures(water height equivalent pressure heads of water and LNAPL) and saturations. In this ﬁgure, the lower curve corresponds to the eﬀective
water saturation distribution in the subsurface, while the upper curve
corresponds to the total liquid (LNAPL and water) eﬀective saturation
distribution.
Porous media with heterogeneous pore systems containing secondary porosity or multiple systems of pores cannot be represented byrelatively standard sigmoidal-type capillary pressure curves. Fig. 1cgives an example of the eﬀective water and total saturation distribution
for a multimodal pore size porous medium. As can be seen, the curves
require adaptations since they are diﬀerent from those having normal

2. Theoretical and experimental methods
First, we investigate whether the scaling of the air-water capillarypressure-saturation relations to obtain air-LNAPL and LNAPL-water
capillary pressure-saturations relations is valid for multimodal pore
distributions. Scaling approaches for air-water capillary pressure-saturation relationships to obtain air-LNAPL and LNAPL-water capillarypressure-saturation relations for LNAPL contamination events were
presented ﬁrst by Parker et al. (1987), and later adapted by API (2007),
among others. Initially, the scaling was used for porous media having
normal or lognormal pore-size distributions. We are not aware of studies where the scaling format has been tested for multimodal
577

Journal of Environmental Management 237 (2019) 576–584

M.A. Alfaro Soto, et al.
k

[

(

) ]

Se[w]m (z ) = ∑ wi 1 + α[nw]i h[nw] Ni −Mi(5)h[nw] = (1 − ρr )(z − z nw )

(6)

k

[

(

) ]

Se[t ] m (z ) = ∑ wi 1 + α[an]i h[an] Ni −Mi
i=1

(7)

σaw
αi
σan

(8)

h[an] = (ρr )(z − z an )(9)α[an]i =or lognormal pore-size distributions. Hence, the LDRM methodology(API, 2007) would not yield accurate predictions of ﬂuid saturations insuch multimodal media since the approach assumes normal or lognormal pore-size distributions (i.e., using parameters from conventional
BC and VG type functions).

where α[an]i is the air-LNAPL scaling factor (API, 2007) and the otherparameters are as deﬁned previously. For elevations below zan,Se[t]m = 1.
Equations (4) and (7) are a function of the multimodal capillarypressure curve, as per the linear overlay of the subcurves adapted from
Durner's model given by Equation (3). The eﬀective water and total
saturations can be represented by the following:

2.3. Adaptations of the eﬀective saturation curves of water and total liquids
Various models can be used for describing capillary pressure-saturation relationships. One frequently used model is the VG equation
given by (van Genuchten, 1980):(1)

where hc is the capillary head [L−1] between the nonwetting andwetting ﬂuids, α [L−1], N [-] and M [-] are VG shape parameters, and Se
is eﬀective wetting ﬂuid saturation deﬁned as

Se [w]m =Sw − Swr
1 − Swr − Snr(10)

Se [t ]m =Sw + Sn − Swr − Snr
1 − Swr − Snr(11)where Sn and Snr correspond to LNAPL saturation and residual LNAPL
saturation, respectively. We note that Equations (10) and (11) are thesame as those used for unimodal eﬀective saturations in API (2007),with Swr and Snr being constants.
The water saturation distribution is obtained by rearranging
Equation (10) to give(2)

Se = (S − Sr )/(1 − Sr )(4)σα[nw]i = aw αi
σnw

in which, as in API (2007), α[nw]i is the LNAPL-water scaling factor, andρr represents the LNAPL speciﬁc gravity (i.e., ρr = ρn/ρw, where ρn andρw are the LNAPL and water mass densities, respectively). For elevations(z) below znw we have full water saturation, and hence Se[w]m = 1.
Similarly, the eﬀective total saturation distribution Se[t]m at static
equilibrium conditions as a function of elevation above zan (the airLNAPL interface elevation in a nearby well) in a multimodal medium
with a continuous LNAPL present, can be written as

Fig. 1. (a) Schematic of the elevation of ﬂuids in a monitoring well (API, 2007).
The elevations zgs, zan, znw and zaw correspond to a datum located at or below
znw, the height of the air-LNAPL interface in the well, the height of the LNAPLwater interface in the well, and the height of the water table without a LNAPL,
respectively, bn represents the thickness of the:NAPL layer in the well, and hn
and hw are the distance between the reference level and zan and zaw, respectively(after Farr et al., 1990); (b) representative distribution of LNAPL in a unimodal
porosity medium; (c) representative distribution of LNAPL in a trimodal porosity medium.

Se = [1 + (αhc ) N ]−M

i=1

where S and Sr correspond to the wetting ﬂuid saturation at a given
capillary head hc and the irreducible (or residual) wetting ﬂuid saturation, respectively. The most common form of equation (1) assumes
that the parameters M and N are related by M = 1-1/N (van Genuchten,1980). API (2006) contains a database listing VG parameters, as well asother physical properties, to help investigators who need to characterize the movement and distribution of a LNAPL at sites where theparameters and physical properties cannot be measured by direct
methods. However, the parameters furnished by API are characteristicof soils with a unimodal pore-size distribution, which makes their use
somewhat limited.
Various authors have presented models for the characterization ofsoils having multimodal pore-size distributions, including Ragab et al.
(1981), Smettem and Kirkby (1990), Othmer et al. (1991), Mallants
et al. (1997), Alfaro Soto et al. (2008), Rudiyanto et al. (2013), Li et al.
(2014). A model similar to equation (1) for soils with a multimodal poresystem was formulated by Durner (1994):

Sw (z )m = Swr + (1 − Swr − Snr ) Se[w]m (z )(12)with Se[w]m given by Equation (4). The LNAPL saturation is obtained ina way similar to that proposed by API (2007). The distribution holds forz < zan, and is given by:

Sn (z )m = 1 − Sw (z ) when z ≤ z an

(13)

In the same way, the LNAPL saturation for z > zan can be obtained
from Equations (10) and (11):Sn (z )m = Snr + (1 − Swr − Snr )(Se [t ] m − Se [w] m ) when z > z an

(14)

Equation (11) indicates that the LNAPL residual saturation (Snr)needs to be known before Sn(z)m can be determined. Following API(2007), the LNAPL residual saturation can be related to the initial
LNAPL saturation via the equation:

k

[

]

Se = ∑ wi 1 + (αi h) Ni −Mi
i=1

Snr (z ) = fr Sni (z )
(3)

(15)

in which fr is referred to as the residual LNAPL factor. This factor relates
residual LNAPL saturation (LNAPL trapped within the formation) with
the initial saturation value at the beginning of the recovery period when
imbibing water starts to displace LNAPL into larger pores. Several
studies such as those by Adamski et al. (2005) and Wadill and Parker(1997) suggest that fr varies between about 0.2 for relatively ﬁne-textured soils to 0.5 for coarse-textured soils.
The LNAPL residual saturation can be obtained from Equations (14)where k represents the number of subsystems making up the total distribution of pores and wi are the weighting factors of the subcurves,which vary from 0 to 1 and sum to unity.
Following Durner (1994), the eﬀective water saturation distribution
Se[w]m for static equilibrium conditions as a function of elevation (z)above znw (the LNAPL-water interface in a nearby well) in a multimodal
pore system with continuous LNAPL present can be written as:578

Journal of Environmental Management 237 (2019) 576–584

M.A. Alfaro Soto, et al.

and (15). If Sni(z) is isolated and z ≤ zan, eﬀective total saturation Se[t]mis equal to 1.0 since at these depths saturation is controlled by the effective water saturation Se[w]m. This yields the following equation:

Sni (h) =(1 − Swr )[1 − Se [w] m (z )]1 − fr Se [w] m (h)when z < z an

(16)

For depths where z > zan, eﬀective total saturation Se[t]m) is lessthan 1, with the initial LNAPL saturation Sni(z) given by:

Sni (h) =(1 − Swr )[Se [t ] m (z ) − Se [w] m (z )]1 − fr [1 + Se [w] m (z ) − Se [t ] m (z )]when z > zan(17)

2.3.1. Free NAPL and recoverable-LNAPL speciﬁc volumes of in theformation
Two important variables when addressing the cleanup of subsurface
LNAPLs are the amount of free LNAPL in the formation (i.e., LNAPL
exceeding residual LNAPL saturation), and how much of this free
LNAPL is retrievable using LNAPL recovery technologies. API (2007)
provides a method for determining free-LNAPL and recoverable-LNAPL
speciﬁc volumes on the basis of the LNAPL elevation in a monitoring
well. In addition to quantifying the generally unknown volumes of
LNAPL, API's approach makes it possible to simulate the ability ofcertain technologies to recover LNAPL releases from the subsurface.
The LDRM software (API, 2007) could possibly be used for this purpose.
Determination of the speciﬁc volume of free LNAPL in soils with amultimodal pore-size distribution requires integration of the area below
the LNAPL saturation curve, considering the total porosity of themedium, i.e.:

Fig. 2. Experimental data of the capillary pressure–saturation curves for airwater (AW) and air-LNAPL (AO) systems (HC: hanging column; FP: ﬁlter
paper).

3.1. Air-water and air-LNAPL capillary pressure curves of multimodal poresystems
Fig. 2 shows capillary pressure-saturation curves obtained for the
air-water and air-LNAPL systems of one of the samples of a sandy loam
soil typical of areas in and near the city of Rio Claro in the state of São
Paulo, Brazil. The air-water data show a small, but noticeable, drop
close to saturation (near a capillary pressure head of about 5 cm) due todesaturation of macropores that are typical of lateritic soils (Alfaro Soto
et al., 2008). The air-water data also show signiﬁcant desaturation at
capillary pressures above about 105 cm. The latter drop, often observed
for tropical soils, is contrary to the S-shaped nature of the traditional
VG function given by Equation (1).
Fig. 3 shows best-ﬁt curves to the air-water data assuming unimodal, bimodal and trimodal pore-size distribution descriptions. The
curves are given by Equation (1), Equation (3) with k = 2, and Equation(3) with k = 3, respectively. Optimizations were carried out using Levenberg-Marquardt, (1963) type nonlinear regression techniques. The
results indicate that the trimodal formulation ﬁts the data better thanthe uni- and bi-modal formulations, especially at the higher capillarypressure heads. The three subcurves of the trimodal model provide
appropriate descriptions of both the initial decrease near 5 cm, theclassical middle part of the curve, and then desaturation at higher
pressure heads. While excellent descriptions are obtained, we acknowledge that this is done at the expense of having more adjustableparameters in the capillary pressure-saturation relationship.

zmax

Dn (bn )m =∫ nSn (z )m dz

(18)

znw

where n represents the porosity of the medium. The recoverable LNAPL
speciﬁc volume (Rn) is given similarly by the area between the LNAPL
saturation curve and the residual saturation curve. This can be expressed as:zmax

Rn (bn )m =∫ n (Sn (z )m − Snr (z ) m) dz

znw

(19)

Evaluation of Equations (18) and (19) requires knowledge of the
maximum elevation (zmax) of free LNAPL as triggered by capillary rise.
zmax may be obtained from Equations (4) and (7), assuming thatSe[w]m (z ) = Se[t ]m (z ) when Sn(z)m = Snr, as follows:k∑ wi {1 + [α[nw]i (1 − ρr )(zmax − znw ) ]Ni }−Mi
i=1
k

{

[

] }

= ∑ wi 1 + α[an]i ρr (zmax − z an ) Ni −Mi
i=1

(20)

Equation (20) can be solved readily for zmax using relatively standard numerical techniques.
3. Results and discussion
We ﬁrst show the measured air-water and air-LNAPL capillarypressure curves of our soils having multimodal pore-size distributions.
Using these data, the ability of the unimodal model used by API (2007)
to reproduce air-LNAPL capillary pressure curves from the air-water
curves is evaluated. Next, we will show results of the new model formulation using bi- and trimodal capillary pressure curves, as well as thecurrent API model (2007), for determining LNAPL speciﬁc volumes forsoils having multimodal pore-size distributions, and discuss the results.

Fig. 3. Experimental and optimized air-water capillary pressure curve data for a
sandy loam sample (HC: Hanging column; FP: Filter paper: VG-U: van
Genuchten unimodal; VG-B: van Genuchten bimodal; VG-T: van Genuchten
trimodal).
579

Journal of Environmental Management 237 (2019) 576–584

M.A. Alfaro Soto, et al.

Fig. 4. Predicted eﬀective saturation of the air-LNAPL system (AO-P) using theunimodal, bimodal and trimodal models (Equations (7)–(9)) versus eﬀective
saturation of the air-LNAPL system (AO) obtained by optimizing the experimental data using the unimodal, and bimodal and trimodal formulations(Equation (3)).

3.2. Comparison of eﬀective saturation of air-LNAPL system obtained
directly and indirectly for multimodal pore systems
The air-water capillary pressure head curves were transformed to
air-LNAPL systems following API (2007) procedures. The air-water andair-LNAPL interfacial tensions were σaw = 71.6 mN m−1 and
σao = 27.6 mN m−1, respectively, and the LNAPL speciﬁc gravity(ρr) = 0.8. Fig. 4 shows a comparison of eﬀective saturations (AO)
obtained by ﬁtting the experimental data using the unimodal (Equation(1)) and multimodal (Equation (3)) models with predicted eﬀective
saturation values (AO-P) using the unimodal, bimodal and trimodal
formulations (the latter two obtained using Equations (7)–(9)). The
data show that the indirect uni- and bi-modal models considerably
overestimate eﬀective saturation values at relatively high capillarypressure heads (eﬀective measured AO saturations between 0.2 and0.5), and slightly overestimate directly measured saturations at relatively low capillary pressure heads. By comparison, the trimodal model
was successful in describing all of the AO eﬀective saturations using thepredicted values.

Fig. 5. Measured capillary pressure–saturation data and optimized curves assuming unimodal, bimodal, and trimodal pore-size distributions. Plots are for:a) Soil a (a sandy loam from Rio Claro, Brazil; b) Soil b (a silty ﬁne sand fromSão Carlos, Brazil); c, d, and e) Soils c, d and e (a sandy loam from Bekkevoort,
Belgium; and f) Soil f (a silty clay from São Carlos, Brazil).
Table 1
Estimated adjustment parameters of the capillary pressure curves.
Parameters

UnimodalSr
α1 (m−1)
N
M
R2
BimodalSr
W1W2
α1 (m−1)α2 (m−1)
N1
N2
M1
M2
R2
Trimodal
Sr
W1
W2
W3
α1 (m−1)α2 (m−1)α3 (m−1)
N1
N2
N3
M1
M2
M3
R2

3.3. Comparison of speciﬁc LNAPL volumes assuming unimodal andmultimodal pore-size distributions
The results in Figs. 2–4 are for a lateritic sandy loam soil from
Brazil. We additionally compared expected subsurface free LNAPL andrecoverable LNAPL distributions as predicted with the unimodal andmultimodal equations assuming hypothetical scenarios of the contaminated aquifers. For this analysis, we used data from six soils(samples a through f) having diﬀerent hydraulic properties. Soil a is asilty sand from Rio Claro (Alfaro Soto et al., 2007), while soils b and fwere clayey sand and silty clay samples, respectively, from São Carlos(Graciosa et al., 2008). These three soils were chosen because of theirmultimodal pore-size distributions typical of many weathered tropicalsoils. In addition, soils c, d and e are undisturbed sandy loam samples
taken from the Bekkevoort ﬁeld site near the city of Leuven, Belgium(Mallants et al., 1997). The latter soils were selected because of theirdiﬀerent hydropedological characteristics, varying from relatively
coarse (a) to very ﬁne-textured (f) media, and the expected diﬀerent
fate of LNAPLs in these media.
Fig. 5 shows the measured air-water capillary pressure data, as wellas the optimized curves based on Equations (1) and (3) for the unimodaland multimodel (bimodal and trimodal) functions, respectively. Valuesof the optimized parameters are listed in Table 1. The results in Fig. 5highlight the diﬃculties of describing the data using only the unimodal
formulation, with the bimodal and especially trimodal models providing better ﬁts to the experimental data. Higher coeﬃcients of

Soil
A

b

c

d

e

f

0.252
2.71
9.68
0.076
0.933

0.291
544.2
9.61
0.0388
0.980

0.077
0.010
0.579
3.72
0.945

0.125
0.690
0.984
0.717
0.971

0.117
0.010
0.498
3.87
0.979

0.224
0.019
0.995
10.07
0.935

0.01250.646
0.354
2.14560.0010
5.71
3.18
0.267
0.807
0.977

0.277
0.152
0.848
0.136
378.2
9.99
4.99
1.000
0.123
0.995

0.021
0.499
0.501
0.1380
11.12
2.94
0.474
1.000
0.409
0.992

0.080
0.920
0.080
1.810
68.4
5.76
2.85
0.0736
0.564
0.978

0.085
0.306
0.694
0.085
3.18
2.91
0.643
1.000
0.465
0.996

0.248
0.499
0.501
0.1275
1.13
9.90
2.95
1.000
0.972
0.999

0.000
0.333
0.097
0.570
0.0010
0.313
1.54
3.03
3.01
3.69
0.608
0.23240.9029
0.977

0.278
0.164
0.201
0.635
0.137
11.5
168. 9
9.99
3.12
4.10
10.0
1.00001.0000
0.998

0.140
0.501
0.332
0.167
0.274
0.078
80.9
1.38
4.44
10.0
1.000
1.000
0.187
0.999

0.170
0.361
0.543
0.096
0.054
0.894
82.0
10.0
4.09
10.0
1.000
1.000
0.1075
1.000

0.000
0.217
0.700
0.083
0.082
3.59
59.0
6.17
1.41
1.000
1.000
0.1511
1.000
0.996

0.250
0.505
0.065
0.430
0.1425
0.152
3.71
3.37
3.15
3.03
1.00000.92030.62630.968

determination (R2) resulted (Table 1) using the trimodal formulation(notably for soils a, b, c, d, and e), with the bimodal formulation also
giving relatively good results for soils a, e, and f. Some of the capillary580

Journal of Environmental Management 237 (2019) 576–584

M.A. Alfaro Soto, et al.

Fig. 7. Calculated LNAPL, water and residual NAPL saturation proﬁles for Soil dconsidering: a) unimodal and b) trimodal pore-size distributions as based on the
ﬁtted optimized capillary pressure-saturation curves.

Fig. 6. Relative pore-size distributions assuming unimodal, bimodal and trimodal pore systems for a) Soil a (a sandy loam from Rio Claro, Brazil; b) Soil b(a silty ﬁne sand from São Carlos, Brazil); c, d, and e) Soils c, d and e (a sandy
loam from Bekkevoort, Belgium), and f) Soil f (a silty clay from São Carlos,
Brazil).

3.4. Predicting water and LNAPL saturation distributions
We now show how the multimodal nature of the pore structure of asoil may aﬀect water and LNAPL saturation distributions in the subsurface. The analysis assumes a hypothetical 4-m deep homogeneous
soil proﬁle, with a water table 2.0 m below the soil surface and amaximum LNAPL thickness of 2.0 m in a nearby well (Fig. 8). For the
ﬂuids we used the following data: LNAPL speciﬁc gravity ρr = 0.8, andinterfacial air-water (σaw), air-LNAPL (σan) and LNAPL-water (σnw)
tensions of 65, 25 and 15 mN m−1, respectively. We also used the optimized air-water hydraulic parameters listed in Table 1. Additionally,
we assumed a total porosity of 0.46 and residual saturation (Snr)varying along the soil proﬁle. For the residual factor (fr), an averagevalue of 0.3 was used, since this coeﬃcient was shown to vary from 0.2
to 0.5 for coarse- and ﬁne-textured soils (API, 2007).
Figs. 7 and 8 show the saturation proﬁles obtained for soil samples dand e, respectively. Each ﬁgure contains a proﬁle obtained using the
API model (2007) for a unimodal pore distribution (Figs. 7a and 8a), as
well as proﬁles based on Equations (12) and (17) for the multimodal
pore distribution (Figs. 7b and 8b). The initial distribution of the saturation curve for soil sample d (Fig. 7b) showed a trimodal distributionof pores, with the saturation curves for water (Sw), LNAPL (Sn) and
residual LNAPL (Snr) being diﬀerent from those obtained for a unimodal
distribution (Fig. 7a). This is not that obvious when comparing thebimodal distribution of Fig. 8b with the unimodal distribution of Fig. 8afor soil sample e. However, for these samples (Figs. 7 and 8), as well as
for samples d, b, c and f, appreciable diﬀerences occur between theunimodal and multimodal saturation curves, even though they are notnecessarily pronounced. This suggest diﬀerences in the magnitude of
the respective LNAPL speciﬁc volumes (free and recoverable), as will beshown later.

pressure curves (notably for soils c, d, and e) revealed rapid changes
near full saturation as a function of capillary pressure, which can belinked to the presence of macropores (Mallants et al., 1997), whileother samples (soils d, e, and b) showed secondary scaling, a possibly
result of aggregation in especially loamy soils (e.g., Smettem and Kirby,1990).
The diﬀerences in the capillary pressure-saturation relations between the samples in Fig. 5 resulted from diﬀerences in their pore-size
distributions. To demonstrate this, we use a simple lognormal distribution of the form (e.g., Durner, 1994):

P (S ) =dS (h)
d log h(21)

where S(h) is described using Equation (1) or (3). The resulting lognormal pore-size distribution for each capillary pressure curve in Fig. 5
is shown in Fig. 6 as a function of pF = log(h), with h in cm. The plots
indicate that the maximum pore-size density for the unimodal distribution ranges from about pF = 0 to mostly pF = 3 (10 m). As wasfound by Mallants et al. (1997), the distributions based on bimodal
capillary pressure curves show, on the left side of the plots, a system of
secondary pores resulting from structuring, probably caused by theaggregation of smaller particles, and on the right, a system of primarypores. The same occurs with a distribution of pores consistent with thetrimodal S(h) or P(h) curves, indicative of a heterogeneous system ofpores. Generally, when a trimodal system is well deﬁned (Fig. 6b,c,d),an initial peak occurs close to saturation reﬂecting macroporosity
contributions.

581

Journal of Environmental Management 237 (2019) 576–584

M.A. Alfaro Soto, et al.

Fig. 10. Recoverable LNAPL speciﬁc volumes of soils a to f for diﬀerent LNAPL
thicknesses.

Many studies (e.g., Schaap and van Genuchten, 2006; Jarvis, 2008;
Araujo et al., 2017, among others) have shown that microporous, aggregated or structured soils generally have much larger hydraulicconductivities at and near full saturation. This implies that the dynamics of water ﬂow will be diﬀerent, and will inﬂuence multiphasic
ﬂow and determination of LNAPL speciﬁc volumes. In an attempt toevaluate possible errors resulting from the use of a unimodal model for
multimodal porosity soils, Dn and Rn were calculated for the diﬀerent
soils, by varying the hypothetical LNAPL thickness (2.0; 1.0, 0.5 m and0.25 m), using the traditional model of API (2007). Associated errors
were calculated using the expression

Fig. 8. Calculated LNAPL, water and residual NAPL saturation proﬁles for Soil econsidering: a) unimodal and b) trimodal pore-size distributions as based on the
ﬁtted optimized capillary pressure-saturation curves.

Error = 100
3.5. Predictions of free LNAPL and recoverable LNAPL speciﬁc volumes for
multimodal pore-size distribution formations

VEn, m − VEn
VEn

(22)

where VEn and VEn,m represent the LNAPL speciﬁc volume (free Dn or
recoverable Rn) calculated for unimodal and multimodal pore distributions, respectively. Fig. 11 shows for each soil the percent error in
Rn as a function of LNAPL thickness. Results indicate that for some ofthe soils Rn will be overestimated, and for others underestimated whena unimodal model is used. Moreover, the magnitude of the error increases with a decrease in the thickness of the LNAPL layer in themonitoring well. This example shows that the calculated Rn values forthe unimodal model underestimated by about 75% and overestimatedby about 200% for the thickest LNAPL layer (2.0 m), and underestimated by about 25% for the thinnest LNAPL layer (0.25 cm).

The free LNAPL (Dn) and recoverable-LNAPL (Rn) speciﬁc volumes
were calculated using Equations (18) and (19), respectively, as well aswith Equation (20) to determine the maximum capillary elevation of
the LNAPL. Figs. 9 and 10 present results of Dn and Rn for the soils
studied, with diﬀerent hypothetical thicknesses of the LNAPL (2.0, 1.0,
0.5 and 0.25 m) in the monitoring well, for the most representative pore
distributions (trimodal for a to e and bimodal for f). The ﬁgures showthat for a given thickness of the LNAPL zone in the monitoring well, thespeciﬁc volumes are not necessarily related only to soil texture. For
example, despite the fact that soil B is not coarse, it has a greater
speciﬁc volume for a given thickness of LNAPL. This is probably a resultof the multimodal nature of its pore-size distribution.

Fig. 9. Free LNAPL speciﬁc volumes of soils a to f for diﬀerent LNAPL thicknesses.

Fig. 11. Errors in calculated recoverable LNAPL speciﬁc volumes (Rn) with theuse of the unimodal model for soils having a multimodal pore-size distribution.
582

Journal of Environmental Management 237 (2019) 576–584

M.A. Alfaro Soto, et al.

4. Conclusions

Durner, W., 1994. Hydraulic conductivity estimation for soils with heterogeneous pore
structure. Water Resour. Res. 30, 211–223.
Farr, A.M., Houghtalen, R.J., McWhorter, D.B., 1990. Volume estimation of light nonaqueous phase liquids in porous media. Gr. Water 28 (1), 48–56.
Gerke, H.H., van Genuchten, M.Th, 1993. A dual-porosity model for simulating the preferential movement of water and solutes in structured porous media. Water Resour.
Res. 29 (2), 305–319.
Graciosa, M.C.P., Mendiondo, E.M., Chaundhry, F.H., 2008. Hydraulic simulation of
rainwater trenches. Revista Brasileira de Recursos Hídricos. 12 (2), 89–99 (inPortuguese).
Hall, R.A., Blake, S.B., Champlin Jr., S.C., 1984. Determination of hydrocarbon thicknesses in sediments using borehole data. In: Proc. 4th Natl. Symp. On Aquifer
Restoration and Ground Water Monitoring. Natl. Ground Water Assoc., Dublin, OH,
pp. 300–304.
Hampton, D.R., Miller, P.D.G., 1988. Laboratory investigation of the relationship between
actual and apparent product thickness in sands. In: Proc. Conf. On Petroleum
Hydrocarbons and Organic Chemicals in Ground Water: Prevention, Detection, andRestoration. Nat. Ground Water Assoc., Dublin, OH, pp. 157–181.
Hernández-Espriú, A., Martínez-Santos, P., Sánchez-León, E., Marín, L.E., 2012. Freeproduct plume distribution and recovery modeling prediction in a diesel-contaminated volcanic aquifer. Phys. Chem. Earth 37–39, 43–51.
ITRC, 2016. Light Nonaqueous - Phase Liquids: Science, Management, and Technology.
Two-Day Workshop, Interstate Technology and Regulatory Council. LNAPLs team,Atlanta, CA 5-6 April.
Jarvis, N.J., 2008. Near-saturated hydraulic properties of macroporous soils. Vadose Zone
J. 7, 1302–1310. https://doi.org/10.2136/vzj2008.0065.
Jeong, J., Charbeneau, R.J., 2014. An analytical model for predicting LNAPL distribution
and recovery from multi-layered soils. J. Contam. Hydrol. 156, 52–61.
Lenhard, R.J., Parker, J.C., 1990. Estimation of free hydrocarbon volume from ﬂuid levels
in monitoring wells. Gr. Water 28 (1), 57–67.
Lenhard, R.J., Oostrom, M., Dane, J.H., 2004. A constitutive model for air-NAPL-waterﬂow in the vadose zone accounting for immobile, non-occluded (residual) NAPL instrongly water-wet porous media. J. Contam. Hydrol. 73, 283–304.
Lenhard, R.J., Rayner, J.L., Davis, G.B., 2017. A practical tool for estimating subsurface
LNAPL distributions and transmissivity using current and historical ﬂuid levels ingroundwater wells: eﬀect of entrapped and residual LNAPL. J. Contam. Hydrol. 205,
1–11.
Lenhard, R.J., Sookhak Lari, K., Rayner, J.L., Davis, G.B., 2018. Evaluating an analyticalmodel to predict subsurface LNAPL distributions and transmissivity from current andhistoric ﬂuid levels in groundwater wells: comparing results to numerical simulations. Ground Water Monit. Remediation 38, 75–84. https://doi.org/10.1111/gwmr.
12254.
Li, X., Li, J.H., Zhang, L.M., 2014. Predicting bimodal soil–water characteristic curves andpermeability functions using physically based parameters. Comput. Geotech. 57,
85–96.
Madi, R., Rooij, G.H., Mielenz, H., Mai, J., 2018. Parametric soil water retention models: a
critical evaluation of expressions of the full moisture range. Hydrol. Earth Syst. Sci.
22, 1193–1219.
Mallants, D., Tseng, P.H., Toride, N., Timmerman, A., Feyen, J., 1997. Evaluation ofmultimodal hydraulic functions in characterizing a heterogeneous ﬁeld soil. J.
Hydrol. 95, 172–199.
Marquardt, D.W., 1963. An algorithm for least-squares estimation of nonlinear parameters. SIAM J. Appl. Math. 11 (2), 431–441.
Mercer, J.W., Cohen, R.M., 1990. A review of immiscible ﬂuids in the subsurface: properties, models, characterization, and remediation. J. Contam. Hydrol. 6, 107–163.
Miguel, M.G., Teixeira, R.S., Padilha, A.C.C., 2006. Soil-water characteristic curves for thelateritic soil of the city of Londrina region. Revista de Ciência e Tecnologia 12 (24),
63–74 (in Portuguese).
Newell, C.J., Acree, S.D., Ross, R.R., Huling, S.G., 1995. Light Nonaqueous Phase Liquids.
U.S. Environmental Protection Agency, Washington, DC EPA/540/S-95/500.
Othmer, H., Diekkruger, B., Kutilek, M., 1991. Bimodal porosity and unsaturated hydraulic conductivity. Soil Sci. 152, 139–149.
Parker, J.C., Lenhard, R.J., 1989. Vertical Integration of Three Phase Flow Equations forAnalysis of Light Hydrocarbon Plume Movement. vol. 5. Transp. Porous Media, pp.
187–206.
Parker, J.C., Lenhard, R.J., Kuppusamy, T., 1987. A parametric model for constitutive
properties governing multiphase ﬂow in porous media. Water Resour. Res. 23,
618–624.
Peters, R.R., Klavetter, E.A., 1988. A continuum model for water movement in an unsaturated fractured rock mass. Water Resour. Res. 24, 416–430.
Ragab, R., Feyen, J., Hillel, D., 1981. Comparative study of numerical and laboratory
methods for determining the hydraulic conductivity function of a sand. Soil Sci. 134,
375–388.
Romano, N., Nasta, P., Severino, G., Hopmans, J.W., 2011. Using bimodal lognormal
functions to describe soil hydraulic properties. Soil Sci. Soc. Am. J. 75, 468–480.
https://doi.org/10.2136/sssaj2010.0084.
Rudiyanto, Toride, N., Sakai, M., van Genuchten, M.Th, 2013. Estimating the unsaturated
hydraulic conductivity of Andisols using the evaporation method. J. Jpn. Soc Phys.
125, 3–15. https://js-soilphysics.com/downloads/pdf/125003.pdf.
Schaap, M.G., van Genuchten, M.Th, 2006. A modiﬁed Mualem-van Genuchten formulation for improved description of the hydraulic conductivity near saturation.
Vadose Zone J. 5, 27–34.
Smettem, K.R.J., Kirkby, C., 1990. Measuring the hydraulic properties of a stable aggregate soil. J. Hydrol. 177, 1–13.
Sookhak Lari, K., Raynera, J.L., Davis, G.B., 2018. Towards characterizing LNAPL remediation endpoints. J. Environ. Manag. 224, 97–105.

The API (2007) model for determining LNAPL speciﬁc and recoverable volumes was modiﬁed for considering soils having a multimodal pore-size distribution. Our modiﬁcations of the API (2007) approach employed a formulation by Durner (1994) for describing
heterogeneous (bi- or tri-modal) pore systems containing both primary(relatively small) and secondary (relatively large) pores. Theoretically,
the approach can describe capillary pressure-saturation relations ofcomplex pore systems involving macropores and/or interaggregate
pores in combination with smaller intra-aggregate and related matrix
pores. We further demonstrated that free LNAPL and recoverable
LNAPL speciﬁc volumes will be diﬀerent depending upon the presenceof unimodal or multimodal (dual-porosity) pore systems. The diﬀerences will likely be greater for cases when the thickness of LNAPL inmonitoring wells is thinner. Our results suggest that estimates of thepermeability to LNAPL, as well as of the performance of certain technologies for recovering free product petroleum liquid released to
groundwater, such as those found in the LDRM software of the API(2007), will depend on the pore structure (i.e., uni- or multi-modal) of
the subsurface. Erroneous predictions of potential LNAPL recoveries via
conventional techniques may result if the correct shapes of the capillarypressure-saturation relations of diﬀerent pore systems are not well deﬁned.
Acknowledgements
The authors would like to thank the Foundation for theDevelopment of UNESP (Fundunesp) for its ﬁnancial support during
research for this paper.
References
Adamski, M., Kremesec, V., Kolhatkar, R., Pearson, C., Rowan, B., 2005. LNAPL in ﬁnegrained soils: conceptualization of saturation, distribution, recovery, and theirmodeling. Ground Water Monit. Remediation 25 (1), 100–112.
AIRE, C.L., 2014. Illustrated Handbook of LNAPL Transport and Fate in the Subsurface.
Contaminated Land: Applications in Real Environments. London, UK. http://www.
claire.co.uk/LNAPL.
Alfaro Soto, M.A., Kumayama, D.M., Chang, H.K., 2007. Calibration of a reﬂectometer forwater ﬂow studies in unsaturated soils. Geociencias 26, 357–368 (in Portuguese).
Alfaro Soto, M.A., Chang, H.K., Vilar, O.M., 2008. Evaluation of fractal scaling of some
Brazilian soils. Rev. Bras. Geociencias 38, 253–262 (in Portuguese).
API, 2006. API Interactive LNAPL Guide. Version 2.0.4. American Petroleum Institute,
Washington, DC.
API, 2007. Light Non-aqueous Phase Liquid Distribution and Recovery Model (LDRM),
Publication 4760. American Petroleum Institute, Washington, DC.
Araújo, A.G., Souza, J.C., Alves, E.C., Silva, I.L.A., Correchel, V., 2017. Application ofmodels of adjustments of soil water characteristics curves in tropical soils. Revista
Mirante, Anápolis (GO) 10 (5), 50–61 (in Portuguese).
ASTM, 2003. Standard Test Methods for Determination of the Soil Water Characteristic
Curve for Desorption Using a Hanging Column, Pressure Extractor, Chilled Mirror
Hygrometer, And/or Centrifuge. Am. Soc. Testing and Materials, Philadelphia, PA
ASTM D6836–02.
ASTM, 2003. Standard Test Method for Measurement of Soil Potential (Suction) Using
Filter Paper. Am. Soc. Testing Materials, Philadelphia, PA ASTM D5298-03.
ASTM, 2014. Standard Test Methods for Surface and Interfacial Tension of Solutions ofPaints, Solvents, Solutions of Surface-Active Agents, and Related Materials. Am. Soc.
Testing and Materials, Philadelphia, PA ASTM D1331–14.
Brooks, R.H., Corey, A.T., 1964. Hydraulic Properties of Porous Media. Hydrol. Colorado
State Univ., Fort Collins, CO Paper No. 3.
Campbell, G.S., Gee, G.W., 1986. Water potential: miscellaneous methods. In: Klute, A.
(Ed.), Methods of Soil Analysis. Part 1, second ed. Agron. Monogr. 9, ASA and SSSA,
Madison, WI, pp. 619–633.
Carducci, C.E., Oliveira, G.C., Zeviani, W.M., Lima, V.M.P., Seraﬁm, M.E., 2013. Bimodal
pore distribution of soils under conservation management system for coﬀee crop.
Eng. Agríc., Jaboticabal. 33 (2), 291–302.
Charbeneau, R.J., 2003. Models for Design of Free-Product Recovery Systems for
Petroleum Hydrocarbon Liquids. API Publication 4729, Am. Petrol. Inst.,
Washington, DC.
Charbeneau, R.J., Johns, R.T., Lake, L.W., McAdams III, M.J., 2000. Free-product recovery of petroleum hydrocarbon liquids. Ground Water Monit. Remediation 20,
147–158.
De Pastrovich, T.L., Baradat, Y., Barthel, R., Chiarelli, A., Fussell, D.R., 1979. Protectionof Ground Water from Oil Pollution. CONCAWE, The Hague, pp. 61.

583

Journal of Environmental Management 237 (2019) 576–584

M.A. Alfaro Soto, et al.

Waddill, D.W., Parker, J.C., 1997. Simulated recovery of light, nonaqueous phase liquidfrom unconﬁned heterogeneous aquifers. Gr. Water 35 (6), 938–947.
Wickramanayake, G.B., Gupta, N., Hinchee, R.E., Nielsen, B.J., 1991. Free petroleumhydrocarbon volume estimates from monitoring well data. J. Environ. Eng. 117,
686–691.
Wijaya, M., Leong, E.C., 2016. Equation for unimodal and bimodal soil-water characteristic curves. Soils Found. 56 (2), 291–300. https://doi.org/10.1016/j. sandf.
2016.02.011.

Spohrer, K., Herrmann, L., Ingwersen, J., Stahr, K., 2006. Applicability of uni- and bimodal retention functions for water ﬂow modeling in a tropical Acrisol. Vadose Zone
J. 5, 48–58.
van Dam, J., 1967. The migration of hydrocarbons in a water bearing stratum. In: Hepple,
P. (Ed.), The Joint Problems of the Oil and Water Industries. Elsevier, Amsterdam, pp.
55–96.
van Genuchten, M.Th, 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.

584