Journal of Contaminant Hydrology 251 (2022) 104089

Contents lists available at ScienceDirect

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

Simulated leaching of PFAS from land-applied municipal biosolids at
agricultural sites
Jeff A.K. Silva a, *, Jennifer L. Guelfo b, Jiří Šimůnek c, John E. McCray daArclight Research & Consulting, Golden, CO, USA
Department of Civil, Environmental, and Construction Engineering, Texas Tech University, Lubbock, TX, United States of AmericacDepartment of Environmental Sciences, University of California, Riverside, United States of AmericadCivil & Environmental Engineering Department, Hydrologic Science & Engineering Program, Colorado School of Mines, Golden, CO, United States of Americab

A R T I C L E I N F O

A B S T R A C T

Keywords:
Land-applied biosolids
PFAS leaching
PFAS vadose zone transport
Numerical modeling
PFAS source term modelingBest management practices

Biosolids are an important resource for agricultural practice but have recently received increased focus as apotential source of per- and polyfluoroalkyl substances (PFAS) in the environment. Few studies have investigated
the transport of PFAS through the unsaturated zone under conditions relevant to biosolids application sites.
Herein, the unsaturated flow and transport model HYDRUS is used to evaluate the leaching of per- and poly­
fluoroalkyl substances (PFAS) from land-applied biosolids used in agricultural practice to determine the impacts
of PFAS leaching on underlying groundwater resources. This numerical case study was based on conditions andoperations at two test sites in central Illinois where biosolids were applied at agronomic rates and where PFAS
contents and desorption characteristics were previously characterized. Each site possessed different vadose zone
soil textural heterogeneity. Simulations were performed under actual present-day meteorological conditions andextended 150 years beyond the initial biosolids application. These long-term simulations demonstrate how soil
equilibrium sorption/desorption processes within the biosolids-amended surface soils effectively control thetransport rate of individual PFAS to groundwater. Air-water interfacial (AWI) adsorption, which is sometimes
considered to be a significant source of PFAS retention in vadose zone soils, was observed to have minimal
impacts on PFAS leaching rates within the biosolids-amended surface soils at these sites. Additionally, the impact
of AWI adsorption was found to be most significant for PFAS transport within the underlying vadose zone soils
when these soils were more texturally homogeneous and considerably less significant within the texturally
heterogeneous soils represented herein. The results of multiple long-term simulations were used to develop an
empirical equation that relates predicted maximum PFAS pore-water concentrations reaching the saturated zone
with changes in PFAS concentrations in the biosolids-amended soil for various biosolids re-application events.
This approach is shown to be very useful in developing site-specific PFAS soil screening levels and/or maximum
leachate levels for PFAS in support of establishing best management practices (BMPs) for land application of
biosolids.

1. Introduction
Biosolids and products derived from biosolids are an important
supplemental source of nutrients for agronomic, silviculture, and land
restoration purposes. Biosolids are used in agriculture because they
contain high concentrations of essential nutrients, primarily nitrogenand phosphorus, which are recycled into crop production to supplement
or replace inorganic fertilizers. Biosolids also replenish soil organic
matter, enhance soil water retention, and improve soil properties tosupport plant root growth. The total mass of biosolids applied to soils in

the United States (U.S.) is substantial; approximately 7.2 million dry
metric tons of biosolids are generated in the U.S. annually, with roughly
55% (or 3.9 million tons) being land applied (Water Environment
Federation WEF, 2017). While the agricultural application of biosolids
represents a beneficial use of an otherwise discarded waste product,
biosolids have also been shown to be a source of heavy metals, phar­maceuticals, and other environmental pollutants (e.g., USEPA, 1990;
Kinney et al., 2006; Harrison et al., 2006; Clarke et al., 2010; Clarke andSmith, 2011; Higgins et al., 2011). This has fueled research aimed atdetermining the fate of these contaminants following land application* Corresponding author.
E-mail address: jsilva@arclightrc.com (J.A.K. Silva).
https://doi.org/10.1016/j.jconhyd.2022.104089
Received 19 July 2022; Received in revised form 12 September 2022; Accepted 23 September 2022
Available online 30 September 2022
0169-7722/© 2022 Elsevier B.V. All rights reserved.

J.A.K. Silva et al.

Journal of Contaminant Hydrology 251 (2022) 104089(e.g., Wilson et al., 1996; Chaney et al., 1996; USEPA, 2000; Yang et al.,2012; Yager et al., 2014) in support of best management practices foragronomic use.
Among the environmental pollutants associated with biosolids are
per- and polyfluorinated substances (PFAS). PFAS are a group of envi­
ronmental pollutants that have become ubiquitous in the environment
due to their persistence and widespread use in consumer and industrial
products and applications (e.g., USEPA, 2009; Ahrens and Bundschuh,
2014; Sharma et al., 2016). Recent research on the occurrence and fateof PFAS in biosolids has demonstrated that these chemicals can bepersistent in soil (e.g., Washington et al., 2010; Sepulvado et al., 2011;
Pepper et al., 2021), can be taken up by plants and edible crops (e.g.,
Higgins et al., 2011; Yoo et al., 2011; Blaine et al., 2013; Wen et al.,2016), and can leach from the soil and potentially impact underlying
groundwater resources (Sepulvado et al., 2011). These findings haveraised recent concerns about environmental and human health impacts
associated with PFAS in land-applied biosolids. Such concerns have led
to regulations of PFAS levels in biosolids in states such as Maine andMassachusetts, which require quarterly monitoring of PFAS in biosolids
destined for land applications (Burns, 2019; Massachusetts Departmentof Environmental Protection (MassDEP), 2022).
Many PFAS are surface-active agents (or surfactants) with the po­tential to accumulate at environmental interfaces in the vadose zone
environment. When released into the vadose zone soils, primary in­terfaces include air-water interfaces (AWI) within the water-unsaturated
soil matrix and soil-water interfaces where PFAS can additionallyassociate with the solid phase (i.e., sorption). Both AWI adsorption andsolid phase sorption combine to slow the rate of PFAS transport. The
significance of PFAS adsorption to AWIs is dependent on the magnitudeof the AWI area (Aaw) within the vadose zone, the surface activity of the
PFAS, and the concentration of PFAS in the solution phase. The
magnitude of Aaw within the vadose zone changes temporally andspatially within the vadose zone environment due to changes in soil
moisture content (θw) as infiltrating rainwater moves vertically down­
ward through the vadose zone. The surface activity of a specific PFAS
molecule varies in response to its structural characteristics, such as the
PFAS alkyl chain length (i.e., carbon number (CN)) and functional group.
Further, the amount of adsorbed PFAS depends on both concentrationand solution ionic strength (Silva et al., 2019). PFAS sorption to modelsands and natural soils has shown to be generally non-ideal (e.g.,moderately non-linear and kinetically or diffusion controlled) (Guelfo
et al., 2020; Brusseau, 2020). These and other complexities associatedwith PFAS vadose zone transport highlight challenges in predicting theenvironmental impacts due to the occurrence of PFAS in land-appliedbiosolids.
Standards of practice for land application of biosolids are based on
the United States Environmental Protection Agency (USEPA) Part 503
Rule (i.e., 40 CFR Part 503, USEPA, 1993), which provides guidance forthe agronomic use of biosolids based primarily on nutrient loading limits
and mitigating the leaching of toxic metals. Currently, Part 503 does notconsider other pollutants like PFAS. Further, the current understandingof PFAS leaching from land-applied biosolids, including the potential toimpact underlying groundwater resources, is limited to field studies, inwhich surficial and depth-discrete soil samples have been collected andanalyzed for PFAS over time (e.g., Washington et al., 2010; Sepulvado
et al., 2011; Pepper et al., 2021). While these studies provide important
information about the impacts of past biosolids application, the datacollected typically represent snapshots of the vertical distribution of
PFAS (i.e., type and mass distribution) resulting from PFAS leaching anddo not provide much information needed to evaluate historical andfuture changes in the vertical distribution of PFAS. Prior field studies
also have not provided data needed to evaluate important causal factors
or environmental impacts, such as the effects of PFAS mass application
rates associated with biosolids application rates, the rate of PFAS
transport within the vadose zone, or the resulting impact on ground­water resources. To perform these evaluations, numerical methods

present a useful tool that can account for the dynamic temporal inter­change of PFAS mass within the vadose zone in response to transient
unsaturated flow conditions resulting from patterns of natural infil­
trating precipitation.
The goal of this work was to evaluate the significance of land-appliedbiosolids as a long-term source of PFAS in agricultural applications
under natural flow conditions, and to assess the potential for PFAS
leaching from these surficial sources to underlying groundwater. Data
collected from two field sites located in central Illinois was used toparameterize a version of the well-known HYDRUS unsaturated flow
and transport model (Šimůnek et al., 2016) that was recently modified tosimulate the transport of PFAS in the vadose zone (Silva et al., 2020).
Here, simulated operational variables include vadose zone hydro­
stratigraphy, biosolids loading rate (BLR), and a number of repeatedannual biosolids applications at a given BLR. Transport-specific vari­ables and conditions were additionally included. As will be demon­
strated, the numerical model was used to evaluate PFAS leaching from asource zone and vertical transport in the underlying vadose zone overextended periods, which in turn allowed for observation of the true
timescales for PFAS transport from the surficial source to the underlying
groundwater table.
2. Summary of site characteristics and model preparation
A detailed description of the methods and data used to prepare andparameterize HYDRUS for unsaturated water flow and PFAS transport isprovided in the Supplemental Information (SI). The SI also includes adetailed discussion of the hydrogeologic characteristics of the test sites
used in this evaluation and biosolids application information pertinent
to this work. The SI also provides details pertaining to the numericalmodeling approach.
Briefly, the two test sites selected for this evaluation are locatedwithin agricultural districts in central Illinois and were a part of a pre­vious study evaluating the impacts of biosolids application on soil
fertility and corn yields (Metropolitan Water Reclamation District ofGreater Chicago (MWRDGC), 2008; Higgins et al., 2011). At these sites,
designated Sites K and W, Class B biosolids were applied to individual
test plots and were incorporated into the topsoil via tilling. Biosolids
loading rates (BLR) varied between individual plots and were applied
over three consecutive years (Section S1.1.4.1).
Notably, the objective of this current work was not to mirror test plot
operations for the purpose of making site-specific predictions ofgroundwater impacts resulting from test operations at this site. Rather,
this study used available data from these sites to parameterize HYDRUS
for general use in evaluating the potential impacts of PFAS leaching togroundwater when biosolids are applied. In some cases, simulations
used BLRs in excess of those actually applied at the test sites. However,simulations were performed to be consistent with actual dates of initial
biosolids application and re-application at these sites. Specific site data
used included historical precipitation and evaporation data (Section
S1.1.1; Fig. S1), site soil types, and hydrostratigraphic data (Section
S1.1.3).
Individual one-dimensional (1D) simulation profiles were con­structed for both test sites (i.e., using HYDRUS) that extended from theground surface to the water table (i.e., ~3 m below ground surface (bgs)
at both sites) and incorporated soil textural changes observed at bothsites. Surficial topsoil and vadose zone soils beneath Site K were ho­mogeneous deposits of silty fine sand of glacial origin (Fig. S2). The
surficial soil at Site W was a loam, and the underlying soils consisted oflayers of silty clay, sandy loam, and medium-grained sands (Fig. S2).
Textural characteristics and unsaturated hydraulic parameters for thesoils are provided in Table S1 and Fig. S3.
Much of the biosolids-specific PFAS data used here was derived fromthe work of Sepulvado et al. (2011), who examined the leaching of PFAS
from soils collected from these same test sites. This includes a linear
biosolids loading rate (BLR, Mg/ha) model used to parameterize PFAS
2

J.A.K. Silva et al.

Journal of Contaminant Hydrology 251 (2022) 104089

mass input in the simulator, i.e.:
SBS,PFAS = m BLR

keeping with the actual initial biosolids application for these sites. The
simulated leaching of PFAS from the applied biosolids is characteristic ofa temporally diminishing source (Fig. 1), with the rate of change in PFAS
concentrations in the leachate (Cw) consistent with the comparative
magnitude of the desorption coefficient (Kdes) assigned to individual
PFAS (Table S3). That is, PFAS with lower CNeach from the biosolids
source faster than those with higher CN. This is consistent with field data
where existing sources show more rapid depletion of lower CN (i.e.,lower molecular weight) PFAS (Sepulvado et al., 2011; Bekele et al.,2020).
The variability in Cw superimposed on PFAS leaching trends
observed for some PFAS demonstrates the response of the solution-phase
concentration to AWI adsorption as infiltrating precipitation changedthe θw and magnitude of the Aaw throughout the simulation period. Notethat increasing θw corresponds with decreaseing Aaw, and vice versa.
This variability is most pronounced for PFNA and PFDA and less pro­nounced for the remaining PFAS. This reflects the variability in thesignificance of AWI adsorption vs. interaction with the soil for a given
PFAS relative to their interaction with the soil. The lower CN PFAS (e.g.,
PFHxA, PFHpA, and PFHxS) are less surface-active, and therefore AWI
adsorption is less significant. PFNA and PFDA are among the mostsurface-active solutes considered here (note the magnitudes of the AWI
adsorption coefficients provided in Table S6). Therefore, in these cases,
AWI adsorption competes more effectively with soil sorption processes
to raise and lower Cw when θw and Aaw change within the PFAS source. In
contrast, while PFDA and PFOS share similar surface activity, the tem­poral trend in Cw for PFOS does not exhibit the same degree of variability
observed for PFDA. PFOS, in this case, is more strongly associated withthe soil phase (i.e., highest Kdes; Table S3) which out-competes AWI
adsorption in the distribution of PFAS mass.
Generally, rates of PFAS desorption from the source zones for the
lower CN PFAS were observed to increase with decreasing BLR over the
15-year simulation period (Figs. S4 and S5). However, for the higher CN
PFASs, the rate of decrease in concentration was much less dependent on
the BLR (Figs. S4 and S5). These results are attributed to differences in
the affinity of the different PFAS for the biosolids-amended soil assorption/desorption processes within the source soils retard the trans­
port of the lower CN PFASs than the higher CN PFASs. Additionally, for agiven BLR, the rate of PFAS leaching was found to be largely insensitive
to the rate constant used (Figs. S4 and S5) even when its value is variedby two orders of magnitude (Table S6). This lack of sensitivity is
attributed to (1) the generally high Kdes of the biosolids-amended soils,(2) the relatively small fraction of rate-limited sites associated withthese same soils (F2 = 0.16 in this case), and (3) the low concentration
gradients associated with low PFAS leachate concentrations, which limit
the driving force for mass transfer to the rate-limited sites. The limited
influence of sorption/desorption kinetics on the leaching of PFAS fromshallow source zones for observed for the 15-year simulations performed
here is consistent with a prior long-term leaching evaluation (Kuntz andGrathwohl, 2009). These results demonstrate that equilibrium sorption
would be an appropriate transport model simplification for Sites K and
W. This result is significant because it suggests that when biosolids areapplied to surface soils at typical agronomic loading rates, predicted
rates of PFAS mass loss from the surficial soil source are more dependent
on an accurate characterization of the equilibrium sorption/desorption
parameters than on non-equilibrium parameters.

(1)

where SBS, PFAS is the biosolids-amended soil concentration (ng/gsoil) andm is the slope of the BLR-SBS, PFAS regression, originally provided by
Sepulvado et al. (2011) (Table S2). Specific PFAS of interest in this work
included five homologous perfluoroalkyl carboxylates (i.e., per­fluorohexanoic acid (PFHxA), perfluoroheptanoic acid (PFHpA), per­fluorooctanoic acid (PFOA), perfluorononanoic acid (PFNA),
perfluorodecanoic acid (PFDA)) and two perfluoroalkyl sulfonates (i.e.,
perfluorohexane sulfonate (PFHxS) and perfluorooctane sulfonate(PFOS)). The perfluorinated tails of these surfactants consist of linear
carbon chains and can be described by their carbon number (CN) (e.g.,
CN for PFHxA/PFHxS = 6 and for PFOA/PFOS = 8).
PFAS leaching from the biosolids was modeled as a two-site kinetic
process that assumed sorption and desorption of PFAS from soil occurs
as both an instantaneous process (i.e., Type 1 sites) and as a slower, ratelimited process (i.e., Type 2 sites) where the fractions of Type I (F1) andType 2 (F2) are used in the model (i.e., F1 + F2 = 1; Section S1.1.4.3).
Desorption rate coefficients (k) used to govern the rate-limited release of
PFAS from agricultural soils were derived from another field study(Washington et al., 2010) (Section S1.1.4.2; Table S3).
Representative soil sorption parameters for the vadose zone soils
underlying the biosolids-amended topsoil were derived from publishedlaboratory experiments performed by Guelfo and Higgins (2013) and
Guelfo et al. (2020). Specifically, PFAS sorption to these vadose zone
soils was represented using either a linear sorption model or a non-linear
Freundlich model, depending on the results of the transport experiments(Section S1.1.4.3; Table S4). Sorption and desorption within vadosezone soils were additionally modeled as a rate-limited process using thesame two-site model used for the biosolids-amended topsoil. Represen­tative F-values and k-values for these soils (Table S5) were derived fromthe work of Guelfo et al. (2020).
AWI adsorption was modeled as a Langmuir-type adsorption process
and the dependence of AWI area on changes in soil moisture content (θw)
was modeled using the Leverett thermodynamic model, also described in
the SI (Section 1.1.4.4) and in the recent literature (Silva et al., 2020,2022). The relationship between Aaw and θw used for the site soils isshown in Fig. S3; AWI adsorption parameters used in these simulations
are provided in Table S6.
Finally, simulations used two monitoring points at both sites. The
first was positioned at the base of the tilled topsoil containing biosolids(i.e., 30 cm bgs), which also represents the base of the PFAS source zone,
to monitor the dynamics of PFAS leaching from the source zone. The
second monitoring point was positioned at the top of the capillary fringe,
representing the first point of contact with the saturated zone. This point
provides information about the potential impacts of PFAS leaching frombiosolids sources on groundwater. Pore-water in the capillary fringe istension-saturated (i.e., hydraulic head is less than atmospheric pressure,or h < 0, where reference pressure is defined as 0 at atmospheric pres­sure). However, because it is fully saturated, significant horizontal flow
first occurs here, which exchanges with flowing groundwater at andbelow the water table (defined as h = 0) (Silliman et al., 2002).
3. Results and discussion
3.1. Source zone dynamics

3.2. Vadose zone hydrodynamics

Generally, the magnitude and rate of loss of PFAS mass from thesource ultimately govern the degree and longevity of the impact onunderlying groundwater resources. Therefore, our first objective was tofocus this evaluation on PFAS leaching from the biosolids applied to thesurface soils. Simulations were performed for a single biosolids appli­
cation event occurring in the Fall of 2004 for Site W and the Spring of2005 for Site K. This timing of the simulated application events is in

Vadose zone hydrodynamics are indicative of the interplay betweenchanges in soil moisture and other soil hydraulic and physical proper­
ties. Infiltrating precipitation increases θw, resulting in a corresponding
increase in unsaturated hydraulic conductivity (Kunsat), the downward
flux of water (q; the volume of water per unit area per time) andcontaminant mass. Changes in θw vary spatially and temporally within
the vadose zone resulting in the transient character of the hydraulic3

J.A.K. Silva et al.

Journal of Contaminant Hydrology 251 (2022) 104089

1000

PFHxA
PFNA
PFOS

1000

PFOA
PFHxS

PFHxA
PFNA
PFOS

100

Cw (ng/L)

Cw (ng/L)

100

PFHpA
PFDA

10

1

PFHpA
PFDA

PFOA
PFHxS

10
1
0.1

0.1

0.01
0

5

10

0

15

5

10

Time (years)

Time (years)

(a)

(b)

15

Fig. 1. Example of simulated PFAS release from both test sites: a) Site K and b) Site W. In both cases, BLR = 5 Mg/ha. Note that year 15 in these plots is the calendar
year 2020.

properties profiles (Fig. 2). Similarly, changes in θw result in concomi­tant changes in Aaw (e.g., Silva et al., 2022), which governs the amountof PFAS mass adsorbed at AWIs, and thereby regulates the mass flux ofcontaminants within the vadose zone.
For Site K, the PFAS source zone and the underlying vadose zone
share the same hydraulic properties but different transport character­istics (i.e., soil sorption/desorption coefficients and rate constants).
Hydraulically, Site K is modeled as a homogeneous system, and the
response observed in the parameters in Fig. 2 reflects this homogeneity.
Kunsat values are periodically very low near the soil surface due to low θw
due to evaporation. However, Kunsat generally ranges between 0.01 and1 cm/day below the surface. This is well below the saturated conduc­tivity (21.2 cm/day) of this silty/sandy soil, and q is thus more variablein this system and more dependent on precipitation variability. This
well-drained soil provides a lower pseudo-steady state moisture content
that yields elevated Aaw values (ranging between 200 and 400 cm2/cm3)
and promotes more significant retention via AWI adsorption in theupper vadose zone relative to Site W (see below).
The vadose zone at Site W is texturally and hydraulically heteroge­
neous. The loamy PFAS source zone retains more soil moisture during
the simulation period, which is shown in the uppermost layer of Fig. 2g.
Generally, θw is shown to increase in time within individual strata,
resulting in independent inter-layer changes in Aaw values. The generalincrease in retained moisture increases Kunsat values within the system,
resulting in a reasonably uniform Kunsat distribution vertically. Due tothe higher overall θw, Kunsat in Site W soils can more readily approach thesaturated hydraulic conductivity. The pseudo-steady state flux of waterin the Site W vadose zone is also generally vertically uniform, although
the water flux is still lower in the Site W vadose zone than in the Site K
vadose zone. However, incidents of greater than average rainfall cansignificantly increase q in this system (i.e., note the significant increasein q at Year 15; Fig. 2f). Given the hydraulic and sorption heterogeneity,
PFAS would be expected to transport more slowly in the vadose zone at
Site W vs. Site K, and the transport rate will be less sharply affected by
precipitation patterns. Leaching of PFAS from the source zone would
also be expected to be slowed due to the hydraulic limitations of theunderlying strata.

simulations, the overall rate of PFAS leaching from the biosolidsamended soil source zone into the vadose zone does not change ateither site (Fig. 3). These results again demonstrate the dominance ofsoil sorption/desorption processes on PFAS leaching from the surficial
source soils in these simulated cases. Note, again, the fluctuating Cw
values in simulations that considered AWI adsorption compared tosmooth changes in Cw values when AWI adsorption was neglected. Thisis because AWI adsorption depends on changes in θw and AWI area,
while soil sorption does not. Infiltrating water increases θw and reduces
AWI area, producing the release of PFAS from AWI to the soil solution.
On the other hand, evaporation decreases θw and increases AWI area,
producing the mass transfer of PFAS from the soil solution to AWI.
The effect of AWI adsorption on the simulated transport of PFAS
within the entire vadose zone is shown in Figs. 4 and 5 for Sites K and W,
respectively. For the less surface-active PFAS, AWI adsorption does not
significantly contribute to transport at either site. However, as the sur­face activity of the solute increases, AWI adsorption plays a more sig­
nificant role in the retention of these solutes. For example, at the morehomogeneous Site K, transport of PFOA is observed to be moderately
accelerated, and the transport of PFOS is significantly accelerated when
AWI adsorption is not considered in the model. The effect is not nearly as
contributive at the texturally heterogeneous Site W. While clear differ­ences exist in the rate of change of PFAS concentrations in the upper partof the vadose zone, after 15 years, PFOA concentrations contacting thesaturated zone are not significantly different whether AWI adsorption isconsidered or not. Likewise, after 15 simulated years, PFOS at Site W has
not contacted the saturated zone and has transported to essentially thesame depth (i.e., ~265 cm) regardless of whether AWI adsorption isconsidered.
Assuming equivalent meteorological inputs, the significance of AWI
adsorption on the overall transport of surface-active solutes is dependent
on three primary factors: 1) the degree of affinity for the solid phase (asdemonstrated for the biosolids amended soils), 2) the surface activity ofthe solute (i.e., an affinity for the AWI), and 3) the magnitude of Aaw as itchanges during infiltration events (i.e., higher Aaw values associatedwith lower water content values provide increased AWI adsorption andretention of surface-active solutes). Factor 3 is strongly influenced by
textural heterogeneity. At Site K, the modeled soil textural homogeneityof the vadose zone soils allows for the establishment of a smaller pseudosteady state θw condition that promotes AWI adsorption in the uppervadose zone (Fig. 2) by maintaining higher Aaw. In the absence of AWI
adsorption, transport occurs much faster because transport is now
dependent only on the hydraulic properties of the media and in­teractions with the solid phase. At Site W, the modeled soil textural

3.3. Significance of PFAS AWI adsorption
AWI adsorption can play a significant role in regulating the rate of
PFAS transport in the vadose zone. However, this was not the casewithin the biosolids-amended source zones of Sites K and W. Regardlessof whether AWI adsorption was or was not considered in these4J.A.K. Silva et al.

5
Journal of Contaminant Hydrology 251 (2022) 104089

Fig. 2. Simulated variability in unsaturated hydraulic conductivity, water flux, moisture content, and air-water interfacial area for Site K (a-d) and Site W (e-h). Fluxes are positive upward and negative downward.

J.A.K. Silva et al.

Journal of Contaminant Hydrology 251 (2022) 104089

1000

1000

Site K, with AWI adsorption

Site K, no AWI adsorption
100

Cw (ng/L)

Cw (ng/L)

100

10

1

PFHxA
PFNA
PFOS

PFHpA
PFDA

10

1

PFOA
PFHxS

0.1

PFHxA
PFNA
PFOS

PFHpA
PFDA

PFOA
PFHxS

0.1
0

5

10

15

0

5

Time (years)

10

15

Time (years)

(b)

(a)
1000

1000

100

100

Cw (ng/L)

Site W, no AWI adsorption

Cw (ng/L)

Site W, with AWI adsorption

10

1

PFHxA
PFNA
PFOS

0.1

PFHpA
PFDA

10

1

PFOA
PFHxS

PFHxA
PFNA
PFOS

0.1
0

5

10

15

0

Time (years)

PFHpA
PFDA

PFOA
PFHxS

5

10

15

Time (years)

(c)

(d)

Fig. 3. Simulated PFAS concentration profiles within the biosolids-amended source zone with and without AWI adsorption at Site K (a-b) and Site W (c-d). Rate
constants were those reported by Washington et al. (2010), as presented in Table S3. BLR = 50 Mg/ha in all cases.

heterogeneity limits the contribution of AWI adsorption because thesoils beneath the biosolids-amended source zone maintain a higher θw
and, as a result, lower Aaw. AWI adsorption is therefore not significant in
the upper Site W vadose zone and does not contribute significantly to theretention of PFAS in this system. At Site W, PFAS transport is dominatedby interactions with the soils and the hydraulics of unsaturated flow.
These results demonstrate how soil textural heterogeneity can, undercertain conditions, play a more important role than AWI adsorption in
determining the magnitude and rate of PFAS transport in the vadosezone, even for more surface-active PFAS (e.g., PFOS in these cases). The
downward coarsening sequence of soil layering present at Site W(Fig. S2) promotes this heterogeneity-imposed hydraulic control on
PFAS transport. Alternatively, if the soils beneath the topsoil werecoarser (e.g., sand), steady state soil moisture would likely be lower
within that underlying coarse layer and Aaw values would be muchhigher (see Aaw(θw) curves, Fig. S3c), and AWI adsorption would play agreater role in PFAS retention. However, hydraulic conditions within the
Site W loam topsoil would still limit the rate and magnitude of leachingin this case because of the intermittent release of PFAS from that source
soil due to the capillary barrier effect (i.e., infiltrating precipitation
would need to build up enough hydraulic head at the loam/sand inter­face to overcome the suction in that layer and initiate flow within theunderlying sand). Further, the flux of water entering the coarser media
would be limited by the overall lower hydraulic conductivity of thesurface soil. Therefore, soil textural heterogeneity should be considered
carefully when simulating PFAS leaching in the vadose zone.

3.4. Leaching of PFAS from a land-applied biosolids source
3.4.1. Single-event biosolids application
To better evaluate PFAS leaching in these systems, PFAS transport in
the vadose zone was simulated for a single biosolids application over alonger simulated time period. An example of the simulated arrival of
PFAS at the top of the capillary fringe at both test sites is presented inFig. 6 for BLR = 50 Mg/ha. Immediately apparent in Fig. 6 are the dif­ferences in the time needed for these PFAS to reach the saturated zone.
The less surface-active PFAS (i.e., PFHxA, PFHpA, PFOA, PFHxS) are
generally less retained by soil interactions and AWI adsorption in thesesystems and therefore move more rapidly through the vadose zone.
However, the higher CN PFAS (i.e., PFNA, PFDA, and PFOS in this case)are more surface-active and sorbing, and are considerably moreretained, requiring many decades to transport to the saturated zone andimpact groundwater resources. The times required to achieve maximum
PFAS concentrations are also remarkable, with PFOS shown to require~60 years to reach a peak concentration at the top of the saturated zone
at Site K and ~109 years at Site W.
Cumulative mass flux (CMF) plots additionally demonstrate the dif­ferences in arrival times for PFAS represented in this work (Fig. 6). CMF
values are presented using hectares as the areal dimension to beconsistent with units used for the BLRs and are useful as input datasets (i.
e., source terms) for groundwater (i.e., saturated zone) modeling efforts.
The time to approach maximal CMF reflects the transport characteristicsof the individual PFAS, and the magnitude of the maximum CMF ach­
ieved reflects the initial PFAS mass loading within the source zone for agiven BLR. For these sites, the maximal cumulative PFAS mass flux
6

J.A.K. Silva et al.

Journal of Contaminant Hydrology 251 (2022) 104089

0

-50

-50

-50

-100

-100

-100

Year 0.5
Year 1
Year 2
Year 3
Year 4
Year 5
Year 10
Year 15

-150

-250

Depth (cm)0Depth (cm)

-150

-200

-250

-250

-300

-300
0

50

-150

-200

100

-300
0

100

200

0

300

100

200

Cw (ng/L

Cw (ng/L)

Cw (ng/L)

(a)

(b)

(c)

PFOA - no AWI adsorption

PFHxS - no AWI adsorption

0

-50

-50

-50

-100

-100

-100

Depth (cm)0Year 0.5
Year 1
Year 2
Year 3
Year 4
Year 5
Year 10
Year 15

-150

-200

-250

-150

0

50

100

-150

-200

-200

-250

-250

-300

-300

300

PFOS - no AWI adsorption0

Depth (cm)Depth (cm)

0

-200

Depth (cm)PFOS - with AWI adsorption

PFOA - with AWI adsorption

PFHxS - with AWI adsorption

-300
0

100

200

300

0

100

200

Cw (ng/L

Cw (ng/L)

Cw (ng/L)

(d)

(e)

(f)

300

Fig. 4. PFAS concentrations in pore-water transported to the groundwater table during the first 15 years since biosolids were applied at Site K (BLR = 50 Mg/ha)
simulated with AWI adsorption (a-c) and without AWI adsorption (d-f). Year 15 in these plots is the calendar year 2020.

entering the saturated zone was 1 to 100 g per hectare of biosolidsamended land but, in some cases, required several decades to reachthat loading.
Expanding this evaluation, simulations were performed for addi­
tional single event biosolids applications at BLRs of 5, 10, 25, 50, and
100 Mg/ha. The results of these simulations for PFOA and PFOS arepresented in Fig. 7. Simulated maximum or “peak” leachate concentra­
tions reaching the saturated zone (Cw, max) for all PFAS evaluated herein(Table S7) increased linearly with increasing BLR for all PFAS evaluated,consistent with the use of a linear loading rate model used to parame­terize PFAS mass input in the simulator (Eq. (1)).
The slope of the BLR-Cw, max regression can be used to directly esti­
mate the maximum BLR for a given site and PFAS. However, because the

PFAS mass associated with the applied biosolids can change from source
to source and over time, it is more practical to represent predicted Cw,
max values as a function of SBS, PFAS (e.g., PFOA and PFOS in Fig. 8). SBS,PFAS values were calculated using Eq. (1), and the values for all PFAS
included in this evaluation are presented in Table S8. The slopes of thelinear regressions of this data are likewise provided in Table S9. These
data can then be used to directly estimate Cw, max values for these sites
and simulated conditions if SBS, PFAS is determined for a given biosolids
application.
3.4.2. Multiple biosolids applications
In practice, biosolids are often applied annually over multi-year cy­
cles, as was the case for test Sites K and W (Higgins et al., 2011).
7

J.A.K. Silva et al.

Journal of Contaminant Hydrology 251 (2022) 104089

0

-50

-50

-50

-100

-100

-100

Year 0.5
Year 1
Year 2
Year 3
Year 4
Year 5
Year 10
Year 15

-150

-250

PFHxS

-300
0

5

10

Depth (cm)0Depth (cm)Depth (cm)

0

-200

-150

-150

-200

-200

-250

-250

PFOA

-300
0

15

Cw (ng/L)

200

0

-50

-50

-50

-100

-100

-100

Depth (cm)0Depth (cm)

0

-250

-150

0

5

10

Cw (ng/L)

15

-150

-200

-200

-250

-250

-300

-300

-300

300

PFOS - no AWI adsorption

0

-200

200

(c)

PFOA - no AWI adsorption

Year 0.5
Year 1
Year 2
Year 3
Year 4
Year 5
Year 10
Year 15100

Cw (ng/L)(b)

PFHxS - no AWI adsorption

-150

PFOS

-300

400

Cw (ng/L)(a)

Depth (cm)PFOS - with AWI adsorption

PFOA - with AWI adsorption

PFHxS - with AWI adsorption

0

100

200

300

0

50

100

Cw (ng/L)

Cw (ng/L)

(e)

(f)

(d)

150

Fig. 5. PFAS concentrations in pore-water transported to the groundwater table during the first 15 years since biosolids were applied at Site W (BLR = 50 Mg/ha)
simulated with AWI adsorption:(a-c) and without AWI adsorption (d-f). Year 15 in these plots is the calendar year 2020.

Therefore, additional simulations were performed to examine the effectof multi-year biosolids applications on maximum PFAS concentrationsin leachate reaching the saturated zone. Lower and more agronomic
BLR’s of 5. 10, and 25 Mg/ha, and annual biosolids re-application pe­
riods (Tapp) of 1, 3, 5, and 10 years were used in these simulations. The
remaining conditions and parameters were the same as in previous
simulations. Examples of the simulated arrival of PFOA and PFOS atdifferent application periods are presented in Fig. 9 for both test sites.
Additionally, simulated Cw, max values for PFOA and PFOS resulting fromthe various application periods and loading rates are provided in Fig. 10.
As shown in Fig. 10, simulated Cw, max was observed to increase
nonlinearly with Tapp, and a 2nd-order polynomial regression equation(forced through the origin) was used to represent the trend in the data.
For the Tapp-Cw, max relationships, the degree of non-linearity increased
with increasing BLR, which is most clearly observed in PFOA results (e.
g., Fig. 10a and c) but is also evident in Fig. 9, where the leading edge of
the PFOA arrival profiles are coincident, which indicates the initial PFAS
arrival times are independent of Tapp for a given BLR.
Similar to increasing the BLR, multiple annual biosolids applications
increase the PFAS Cw, max reaching the groundwater table. However,multiple annual applications ultimately provide a greater Cw, max values
for a given BLR due to the re-application of a fresh source of PFAS that
increases the PFAS mass within the surface soils (i.e., SBS, PFAS) and in theleachate. These results indicate a practical limit for biosolids
8J.A.K. Silva et al.

Journal of Contaminant Hydrology 251 (2022) 104089

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

PFHxA
PFNA
PFOS

0

PFHpA
PFDA

Cumulative Mass Flux (mg/ha)

Cw/C w,0Site K
PFOA
PFHxS

1000000
100000
10000
1000
100
10
1
0.1
0.01
0.001
0.0001

0

PFHpA
PFDA

Cumulative Mass Flux (mg/ha)

Cw/C w,0

(b)

(a)
PFHxA
PFNA
PFOS

PFOA
PFHxS

Time (years)

Time (years)

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

PFHpA
PFDA

0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150

10 20 30 40 50 60 70 80 90 100 110 120 130 140 150

Site W

PFHxA
PFNA
PFOS

PFOA
PFHxS

1000000
100000
10000
1000
100
10
1
0.1
0.01
0.001
0.0001

PFHxA
PFNA
PFOS

PFHpA
PFDA

PFOA
PFHxS

0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150

10 20 30 40 50 60 70 80 90 100 110 120 130 140 150

Time (years)

Time (years)

(c)

(d)

Fig. 6. Simulated PFAS transport within the Site K and Site W vadose zones for a single biosolids application with BLR = 50 Mg/ha: (a/c) PFAS relative concen­
trations at the capillary fringe monitoring location, and (b/d) cumulative mass flux of PFAS impacting saturated zone.

applications in agriculture would need to be derived to be protective ofunderlying groundwater resources.
Like regression equations presented in Figs. 7 and 8, those provided
in Fig. 10 could also be used to estimate Cw, max for PFAS at different
values of Tapp and for the BLR and associated SBS, PFAS values applied.
However, a unifying predictive relationship is preferred to include theimpacts of both SBS, PFAS and Tapp for an individual site. To this end, theslope values of the linear SBS, PFAS-Cw, max regressions (Table S9) wereregressed against Tapp and fit using a 2nd-order polynomial regression
equation (forced through the origin). These m- Tapp relationships werethen used in combination with SBS, PFAS-Cw, max regression equations toderive the following polynomial expression that allows the prediction ofCw, max for various combinations of SBS, PFAS and Tapp as:
( (
)2
(
))
Cw,max = SBS,PFAS a1 Tapp + a2 Tapp(2)

3.5. Implications for groundwater contamination
The USEPA recommended the use of 70 ng/L for PFOA and PFOS (i.
e., individual or combined) as the preliminary remediation goal (PRG)
for groundwater contaminated by PFAS, which is a current or potential
source of drinking water for sites where state or local metrics have notbeen adopted (USEPA, 2019). If this same metric is used to definemaximum allowable leaching levels for PFOA and PFOS, the allowableleaching level is shown to have been exceeded at Site K and Site W for asingle BLR of 50 Mg/ha (Fig. 8). Likewise, this PRG is met or exceeded
for several multi-year applications shown in Figs. 9 and 10.
However, more typically, the PRG of 70 ng/L would be for PFOA andPFOS arriving or produced at a down-gradient receptor point (e.g., asurface water body or groundwater well). Therefore, depending on the
distance between the receptor and the biosolids application area,
attenuation processes may allow PFAS leachate concentrations imme­
diately beneath the biosolids application site to be in excess of 70 ng/L.
Well-known analytical solutions for solute transport (e.g., Domenico,1987) have been incorporated into various analytical models to allow
the estimation of dilution attenuation factors (DAF) at user-defined
distances from the source of contamination (e.g., Nevin et al., 1997;
USEPA, 2000). Using the Domenico model for a sorbing and nondecaying solute (Section S3) and relevant site-specific inputs for thesesites, DAF values were calculated for PFOA and PFOS at Site K (Table 2).
The DAF value was used to determine Cw, max for a 70 ng/L arrival at agiven receptor distance from the source. In turn, these Cw, max values
were used to determine the corresponding SBS, PFAS values for PFOA andwhere a1 and a2 are regression coefficients for the PFAS of interest(Table 1). This approach was used for both sites. As proposed, Eq. (2) issite-specific and additionally dependent on the linearity of SBS, PFAS withchanges in BLR for the biosolids applied. Therefore, additional model
simulations would need to be performed for a given site to derive newsite-specific a1 and a2 parameters. However, once derived, Eq. (2) could
be applied to the development of site-specific PFAS soil screening levels
and/or maximum leachate levels for PFAS for proposed SBS, PFAS and Tappto develop best management practices (BMPs) for land application of
biosolids.

9J.A.K. Silva et al.

Journal of Contaminant Hydrology 251 (2022) 104089

Site K300

PFOA - BLR25PFOA - BLR50PFOA - BLR100

250

250

PFOS - BLR25
PFOS - BLR50PFOS - BLR100

200

Cw,max (ng/L)

Cw (ng/L)

200
150
100

y = 1.9267x
R² = 1

150
100
50

50

0

0
0

50

Time (years)

300

0

150

20

40

60

80

100

120

BLR (Mg/ha)(b)

PFOA - BLR25PFOA - BLR50PFOA - BLR100

250

100

(a)

Site W

180

PFOS - BLR25
PFOS - BLR50PFOS - BLR100

y = 1.6018x
R² = 0.9997

PFOA
PFOS

160
140

Cw,max (ng/L)

200

Cw (ng/L)

y = 2.2462x
R² = 1

PFOA
PFOS

150
100

120

y = 1.5208x
R² = 0.9999

100
80
60
40

50

20
0

0
0

50

Time (years)

100

0

150

20

40

60

80

100

120

BLR (Mg/ha)

(d)

(c)

Fig. 7. Results of varied BLRs for a single application event: a/c) Simulated PFOA and PFOS arrival profiles within the Site K and Site W vadose zones, respectively,
and b/d) maximum PFOA and PFOS concentrations contacting the saturated zone for Site K and Site W, respectively.

180

250
y = 128.45x
R² = 1140

y = 11.344x
R² = 1150

Cw,max (ng/L)

Cw,max (ng/L)

200

100

y = 8.09x
R² = 0.9997

120
100
80
60
40

50
0

y = 101.39x
R² = 0.9999

160

PFOA
PFOS
0

5

10

15

20

PFOA
PFOS

20
0
25

0

5

10

15

20

25

SBS,PFAS (μg/kg)

SBS,PFAS (μg/kg)

(b)

(a)

Fig. 8. Maximum PFOA and PFOS concentrations contacting the saturated zone at both sites: a) Site K and b) Site W.

PFOS using the linear regression and slope values described above(Table S9). These SBS, PFAS values would then be the maximum values, or
potential site-specific soil screening levels (SSL) for PFOA and PFOS,
that represent the maximum PFAS mass per unit mass of soil that could

be applied to the land surface either in a single event or accumulated orover multiple application events to avoid exceeding the 70 ng/L PRG at a
distant receptor. For example, a receptor positioned 500 ft (ft) down­
gradient of Site K would see PFOA and PFOS concentrations exceed the10

J.A.K. Silva et al.

Journal of Contaminant Hydrology 251 (2022) 104089

Site K – BLR = 5 Mg/ha
160

PFOA - 1YRPFOA - 3YR
PFOA - 5YR
PFOA - 10YR

140
120

160

PFOS - 1YR
PFOS - 3YR
PFOS - 5YR
PFOS - 10YRPFOA - 1YRPFOA - 3YR
PFOA - 5YR
PFOA - 10YR

140
120

100

Cw (ng/L)

Cw (ng/L)

Site W – BLR = 5 Mg/ha

80
60

100
80
60

40

40

20

20

0

PFOS - 1YR
PFOS - 3YR
PFOS - 5YR
PFOS - 10YR

00

10 20 30 40 50 60 70 80 90 100 110 120 130 140 150

010 20 30 40 50 60 70 80 90 100 110 120 130 140 150

Time (years)

Time (years)

(a)

(b)

Fig. 9. Simulated PFOA and PFOS arrival profiles for different annual application periods at Site K(a) and Site W (b) for BLR = 5 Mg/ha.

700

350
BLR = 5 Mg/ha
BLR = 10 Mg/ha
BLR = 25 Mg/ha

250

Site K - PFOA

200

y = -0.7687x2 + 20.66x
R² = 0.9999

150
100
50
0

2

4

6

8

10

y = 0.2061x2 + 55.84x
R² = 1500

Site K - PFOS

400

y = 0.1186x2 + 21.668x
R² = 1

300
200
100

y = -0.3796x2 + 10.255x
R² = 0.9999

0

BLR = 5 Mg/ha
BLR = 10 Mg/ha
BLR = 25 Mg/ha

600

Cw,max (ng/L)

Cw,max (ng/L)

300

y = -1.8858x2 + 51.363x
R² = 0.9999

y = 0.062x2 + 10.697x
R² = 10120

2Number of Annual Applications

4

(a)

8

10

12

(b)

350

400
BLR = 5 Mg/ha
BLR = 10 Mg/ha
BLR = 25 Mg/ha

300
250

300

Site W - PFOA

200
y = -0.6574x2 + 18.064x
R² = 0.9894

150

BLR = 5 Mg/ha
BLR = 10 Mg/ha
BLR = 25 Mg/ha

350

y = -1.0609x2 + 39.212x
R² = 1

Cw,max (ng/L)

Cw,max (ng/L)

6

Number of Annual Applications

100

y = -0.2098x2 + 38.598x
R² = 1Site W - PFOS

250

y = -0.0591x2 + 14.791x
R² = 1

200
150
100

50

50

y = -0.3218x2 + 8.9217x
R² = 0.9896

0
0

2

4

6

8

10

y = -0.0257x2 + 7.1874x
R² = 10
120Number of Annual Applications

2

4

6

8

10

12

Number of Annual Applications

(c)

(d)

Fig. 10. Dependence of Cw, max values for PFOA and PFOS on the number of annual biosolids applications at Site K (a/b) and Site W (c/d).

70 ng/L limit when SBS, PFOA and SBS, PFOS at Site K exceed 4.32 and 48.9μg/kg, respectively, or when the cumulative BLR in each case exceeds
288 and 247 Mg/ha, respectively. Given the significant amount of timeneeded for PFOA and PFOS transport in the vadose zone and the arrivalof the peak concentration at the saturated zone, the potential for PRG
exceedances at potential receptors should additionally be considered inboth the nearer term (i.e., based on PFOA peak arrival times) and in themuch longer-term (i.e., based on peak PFOS peak arrival), as maximum
solution concentrations at the receptors could occur many decades after
initial biosolids application.
Regarding the continued and long-term use of biosolids at these sites,
the simulated results of this work suggest biosolids application should bea strictly managed activity in a manner consistent with the management
methods outlined in USEPA 40 CFR Part 503 for metals. The BLR values
11

J.A.K. Silva et al.

Journal of Contaminant Hydrology 251 (2022) 104089

Table 1Results of multi-year biosolids application simulations and regression parameters for PFOA and PFOS used in developing and applying Eq. (2).
Site K

Cw,max (ng/L) for different BLR

SSB,PFAS-Cw,max regression

a2

R2

1
1
1
1

− 5.04

137

1

1
1
1
1

0.5915148.1

1Tapp (years)

BLR = 5 Mg/ha

BLR = 10 Mg/ha

BLR = 25 Mg/ha

m

R

PFOA

1
3
5
10

9.41
27.1
42.2
64.5

19.0
54.6
84.8
130

47.9
138
210
325

127.4
366.8
560
866.1

PFOS

1
3
5
10

10.8
32.9
54.8
113

21.8
66.5
111
229

55.5
170
284
579

147.6
450.6
755
1540

Site W

Cw,max (ng/L) for different BLR

m-Tapp regression
a1

PFAS

2

SSB,PFAS-Cw,max regression

m-Tapp regression

PFAS

Tapp (years)

BLR = 5 Mg/ha

BLR = 10 Mg/ha

BLR = 25 Mg/ha

m

R2

a1

a2

R2

PFOA

1
3
5
10

7.24
21.1
39.4
56.6

14.6
42.7
79.6
114

37.2
109
170
286

98.8
288
465
762

1
1
1
1

− 3.09

107.2

1

PFOS

1
3
5
10

7.21
21.3
35.3
69.3

14.8
43.8
72.5
142

38.4
114
188
365

102
302
497
968

1
1
1
1

− 0.53102.1

1

Table 2
DAF and soil screening level calculations for PFOA and PFOS for site K.
Distance from source (ft)DAF

Cw,max

SSB,PFOA

SSB,PFOSPFOA BLR (Mg/ha)

PFOS BLR (Mg/ha)

1
10
100
250
500
750
1000
5280

1
1
1.46
3.25
7.92
15.6
26.4
691

70
70
102
228
554
1095
1848
48,370

0.54
0.54
0.80
1.77
4.32
8.52
14.4
377

6.17
6.17
9.01
20.1
48.9
96.5
163
4265

36
36
53
118
288
568
959
25,105

31
31
45
101
247
487
823
21,534

Note: Cw,max are those calculated for 70 ng/L concentration at a given receptor distance from the source.

greater than those in Table 2 would result in exceedances of the 70 ng/Lleachate level for PFAS at different receptor distances at Site K. How­ever, given the longevity of biosolids as a source of PFAS and the sig­
nificant times involved in PFAS contamination of groundwater for thesesites, the same BLR/SBS, PFAS values may also be the practical upper
lifetime limit of biosolids application for these sites. Concentrations of
PFAS of current regulatory concern (e.g., PFOA and PFOS) in municipal
biosolids will vary between source municipalities and could ultimately
diminish over time as these compounds are phased out of consumer
products. Additionally, methods to treat the biosolids to reduce PFAS
concentrations prior to land application should be evaluated. Any
reduction of the PFAS mass loading will extend the beneficial and safeuse of these products for agricultural applications.
The USEPA recently updated its lifetime drinking water healthadvisory limit for PFOA and PFOS to 0.004 ng/L and 0.02 ng/L,
respectively (USEPA, 2022). However, at the time of this publication, it
is unclear whether the PRG has also been updated, to what degree theserecommended advisory levels will be adapted as state-level regulatory
limits, and whether analytical methods will be developed that canroutinely and reliably evaluate these levels in drinking water. Regard­
less, the approach and methods described in this current evaluation
could prove useful to agricultural practitioners/management as theseenvironmental regulatory metrics change in the future to reflect new
toxicity data and regulation of additional PFAS.

4. Conclusions
Because PFAS are present in land-applied municipal biosolids, it isimportant to understand the fate and transport of these contaminants inthe vadose zone to assess the potential for groundwater contamination.
In this work, numerical simulations were employed to evaluate thesignificance of land-applied biosolids as a long-term source of PFAS in
agricultural applications under natural infiltration conditions using flow
and PFAS-specific transport parameters derived from existing site in­formation or the literature. Simulations were performed for two testsites located in the same climatic zone with differing vadose zonehydrostratigraphy. Results of these simulations demonstrated thefollowing:• Within the biosolids-amended PFAS source, soil sorption/desorption
processes were the dominant retention mechanisms for PFAS, with
rate-limited desorption contributing minimally to the overall release
rate of PFAS from the source. This result is significant because it
suggests that over the long-term, rates of PFAS mass loss from thesurficial biosolids-amended soil source are most dependent on theaccurate characterization of the equilibrium soil interaction param­
eters than on those for non-equilibrium processes. This, in turn,simplifies leaching calculations.
• While AWI adsorption can be an important process for retaining
PFAS in the vadose zone, this process did not contribute significantly
to the overall release rate of PFAS from the biosolids-amended source
soils at either test site. Climate conditions at these sites provided
12

J.A.K. Silva et al.

Journal of Contaminant Hydrology 251 (2022) 104089

persistently elevated soil moisture conditions and correspondingly
lower Aaw values within these surface soils that minimized thecontribution of AWI adsorption.
• Within the vadose zone soils underlying biosolids-amended source
soils, AWI adsorption was observed to be most contributive to theretention of PFAS at Site K, where the vadose zone was characterized
and modeled as a predominantly homogeneous vadose zone of silty
sand. At Site K, lower moisture conditions, and thus larger Aaw
values, were maintained throughout the 150-year simulation period.
This condition allowed AWI adsorption to play a more significant
role in limiting the rate of PFAS transport to groundwater, particu­
larly for the more surface-active PFAS.
• In contrast, textural heterogeneities and finer-grained soils at Site W
collectively minimized the role of AWI adsorption, resulting in fastertransport to groundwater for the same PFAS source condition, evenwith an overall lower unsaturated hydraulic conductivity. Finergrained soils comprising the upper vadose zone at Site W retained
more moisture throughout the period of simulation, resulting inlower Aaw values and less overall retention of PFAS. These results are
consistent with prior work highlighting the impact of vadose zonehydrostratigraphy on the significance of AWI adsorption as amechanism for PFAS retention during leaching from a surficial
source (Silva et al., 2020).
• Individual PFAS leach from a biosolids source, transport through thevadose zone at different rates, and reach the saturated zone atdifferent times. For linear perfluoroalkyl PFAS employed herein, therate of transport is dependent on the PFAS CN, which impacts its
surface activity and the interaction with the soil solid phase and its
surface activity. For two of the PFAS of federal regulatory concern
used here (i.e., PFOA and PFOS), peak concentrations arrive at thesaturated zone many decades apart (Fig. 7). This requires treating
PFOA and PFOS separately when applying USEPA guidance on PFAS
limits for groundwater.
• Increasing the BLR at both sites increased PFAS concentrations in theleachate but did not impact the rate of release from the source zone.
Increasing the number of annual biosolids application events also
increased PFAS concentrations but did not significantly affect therate of transport of individual PFAS through the vadose zone.
However, it did increase the overall mass of PFAS leaching throughthe soil for a given BLR.
• Peak PFAS concentrations reaching the saturated zone increased
linearly with increasing BLR and increased nonlinearly withincreasing biosolids re-application at a given BLR. The results ofthese simulations were used to develop an empirical equation (Eq.
(2)) that was, in turn, used to directly evaluate the impact of
biosolids-amended soil PFAS loading rates and annual re-application
periods on PFAS leaching concentrations for a given test site. Thisapproach could be useful in developing site-specific soil screening
levels for PFAS in support of BMPs for land application of biosolids
containing PFAS at additional sites. Likewise, an up-scaled version ofthis approach could be used to develop state- or county-level soil
screening levels for generic soil conditions most common in areas
where biosolids land application sites are present.
• Example soil screening calculations for Site K (Table 2) also
demonstrate that the numerical modeling approach used here can
support the continued safe use of biosolids with respect to managing
impacts related to PFAS and groundwater contamination in a mannerconsistent with the methods outlined in USEPA 40 CFR Part 503 forother biosolid-specific contaminants.

et al., 2021). For the current case, this additional source of PFAS could
impact the magnitude of Cw, max values and breakthrough profiles in
Figs. 7 and 10 would plateau at a different Cw, max for those PFAS being
produced. Given that the rates and duration of PFAS production can becharacterized, such a source of PFAS could be simulated as an in situ
production process within the source zone, as was done in our priorwork (Silva et al., 2020). However, simulated PFAS arrival times and thegeneral approach proposed here to estabilish soil-screening levels andassist management of PFAS loading rates for biosolids application
remain appropriate.
CRediT 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. Jennifer L. Guelfo: Methodology, Formal analysis,
Resources, Writing – review & editing. Jiří Šimůnek: Methodology,
Software, 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 influence
the work reported in this paper.
Data availability
The data that has been used is confidential.
Acknowledgments
PFAS-specific modifications of the HYDRUS simulator used in thiswork were supported by The US Army Corps of Engineers and The
Department of Defense Strategic Environmental Research and Devel­
opment Program under Contract No. W912HQ21C0040.
Appendix A. Supplementary data
Supplementary data to this article can be found online at https://doi.
org/10.1016/j.jconhyd.2022.104089.
References
Ahrens, L., Bundschuh, M., 2014. Fate and effects of poly- and perfluoroalkyl substancesin the aquatic environment: a review. Environ. Toxicol. Chem. 33 (9), 1921–1929.
Bekele, D.N., Liu, Y., Donaghey, M., Umeh, A., Arachchige, C.S.V., Chadalavada, S.,
Naidu, R., 2020. Separation and lithological mapping of PFAS mixtures in the vadosezone at a contaminated site. Front. Water 2, 597810.
Blaine, A.C., Rich, C.D., Hundal, L.S., Lau, C., Mills, M.A., Harris, K.M., Higgins, C.P.,
2013. Uptake of perfluoroalkyl acids into edible crops via land applied biosolids:
field and greenhouse studies. Environ. Sci. Technol. 47 (24), 14062–14069.
Brusseau, M.L., 2020. Simulating PFAS transport influenced by rate-limited multiprocess retention. Water Res. 168, 115179.
Burns, D., 2019. Memorandum to Licensed Facilities that Land Apply, Compost, orProcess Sludge in Maine, Re: Requirement to Analyze for PFAS Compounds. Maine
Department of Environmental Protection.
Chaney, R.L., Ryan, J.A., O’Connor, G.A., 1996. Organic contaminants in municipal
biosolids: risk assessment, quantitative pathways analysis, and current research
priorities. Sci. Total Environ. 185, 187–216.
Clarke, B.O., Smith, S.R., 2011. Review of ‘emerging’ organic contaminants in biosolids
and assessment of international research priorities for the agricultural use of
biosolids. Environ. Int. 37, 226–247.
Clarke, B.O., Porter, N.A., Marriott, P.J., Blackbeard, J.R., 2010. Investigating the levels
and trends of organochlorine pesticides and polychlorinated biphenyl in sewage
sludge. Environ. Ein. 36, 323–329.
Domenico, P.A., 1987. An analytical model for multidimensional transport of a decaying
contaminant species. J. Hydrol. 91, 49–58.
Guelfo, Jennifer L., Higgins, Christopher P., 2013. Subsurface transport potential of
perfluoroalkyl acids at aqueous film-forming foam (AFFF)-impacted sites. Environ.
Sci. Technol. 47 (9), 4164–4171.

In this work, it was assumed that a finite mass of PFAS was applied aseither a single application or as incremental additions associated with anannual biosolids application. However, recent investigations haveshown the potential for degradation of PFAS precursor compounds
within land-applied biosolids that could be represented as an additionalsource of some PFAS over the time-scales represented herein (Roehler
13

J.A.K. Silva et al.

Journal of Contaminant Hydrology 251 (2022) 104089
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.
USEPA, 1990. National sewage sludge survey; availability of information and data, andanticipated impacts of proposed rule of proposed regulations, proposed rule. Fed.
Regist. 55 (218), 47210–47283.
USEPA, 1993. Standards for the Use or Disposal of Sewage Sludge (40 Code of Federal
Regulations Part 503). Washington, D.C.
USEPA, 2000. Biochlor: Natural Attenuation Decision Support System User’s Manual,
Version 1.0, EPA/600/R-00/008.
USEPA, 2009. Long-Chain Perfluorinated Chemicals (PFCs) Action Plan. Retrieved
October 2018 from. https://www.epa.gov/assessing-and-managing-chemicals-under-tsca/long-chain-perfluorinated-chemicals-pfcs-action-plan.
USEPA, 2019. Interim Recommendations for Addressing Groundwater Contaminatedwith PFOA and PFOS. Fact Sheet., December 2019. https://www.epa.gov/sites/defa
ult/files/2019-12/documents/pfas_groundwater_fact_sheet.pdf (accessed: 14
September 2021).
USEPA, 2022. Technical Fact Sheet: Drinking Water Health Advisories for Four PFAS(PFOA, PFOS, GenX Chemicals, and PFBS). June 2022. https://www.epa.gov/syste
m/files/documents/2022-06/technical-factsheet-four-PFAS.pdf.
Washington, J.W., Yoo, H., Ellington, J.J., Jenkins, T.M., Libelo, E.L., 2010.
Concentrations, distribution, and persistence of perfluoroalkylates in sludge-applied
soils near Decatur, Alabama, USA. Environ. Sci. Technol. 44 (22), 8390–8396.
Water Environment Federation (WEF), 2017. Q&A/Fact Sheet: Land Application of
Biosolids: Human Health Risk Assessment Related to Microconstituents. WSEC-2017FS-014, August 2017. https://www.wef.org/globalassets/assets-wef/3-resources/to
pics/a-n/biosolids/technical-resources/wef-fact-sheet-microconstituents-v25-aug2017.pdf.
Wen, B., Wu, Y., Zhang, H., Liu, Y., Hu, X., Huang, H., Zhang, S., 2016. The roles ofprotein and lipid in the accumulation and distribution of perfluorooctane sulfonate(PFOS) and perfluorooctanoate (PFOA) in plants grown in biosolids-amended soils.
Environ. Pollut. 216, 682–688.
Wilson, S.C., Duarte-Davidson, R., Jones, K.C., 1996. Screening the environmental fate oforganic contaminants in sewage sludges applied to agricultural soils: 1. The potential
for downward movement to groundwaters. Sci. Total Environ. 185, 45–57.
Yager, T.J.B., Furlong, E.T., Kolpin, D.W., Kinney, C.A., Zaugg, S.D., Burkhardt, M.R.,
2014. Dissipation of contaminants of emerging concern in biosolids applied to nonirrigated farmland in eastern Colorado. J. Am. Water Resour. Assoc. 50 (2), 343–357.
https://doi.org/10.1111/jawr.12163.
Yang, Y.-Y., Gray, J.L., Furlong, E.T., Davis, J.G., ReVello, R.C., Borch, T., 2012. Steroid
hormone runoff from agricultural test plots applied with municipal biosollids.
Environ. Sci. Technol. 46 (5), 2746–2754. https://doi.org/10.1021/es203896t.
Yoo, H., Washington, J.W., Jenkins, T.M., Ellington, J.J., 2011. Quantitative
determination of perfluorochemicals and fluorotelomer alcohols in plants frombiosolid-amended fields using LC/MS/MS and GC/MS. Environ. Sci. Technol. 45(19), 7985–7990.

Guelfo, J.L., Wunsch, A., McCray, J., Stults, J.F., Higgins, C.P., 2020. Subsurface
transport potential of perfluoroalkyl acids (PFAAs): column experiments andmodeling. J. Contam. Hydrol. 233, 103661.
Harrison, E.Z., Oaks, S.R., Hysell, M., Hay, A., 2006. Organic chemicals in sewage sludge.
Sci. Total Environ. 367, 481–497.
Higgins, C.P., Paesani, Z.J., Chalew, T.E.A., Halden, R.U., Hundal, L.S., 2011. Persistenceof triclocarban and triclosan in soils after land application of biosolids andbioaccumulation in Eisenia foetida. Environ. Toxicol. Chem. 30 (3), 556–563.
Kinney, C.A., Furlong, E.T., Zaugg, S.D., Burkhardt, M.R., Werner, S.L., Cahill, J.D.,
Jorgensen, G.R., 2006. Survey of organic wastewater contaminants in biosolids
destined for land application. Environ. Sci. Technol. 40 (23), 7207–7215. https://doi.org/10.1021/es0603406.
Kuntz, D., Grathwohl, P., 2009. Comparison of steady-stae and transient fow condition
on reactive transport of contaminants in the vadose zone. J. Hydrol. 369, 225–233.
Massachusetts Department of Environmental Protection (MassDEP), 2022. Information
on MassDEP Efforts to Establish Standards for Per- and Polyfluoroalkyl Substances(PFAS) in Residuals. https://www.mass.gov/info-details/pfas-in-residuals (accessed2022-03-19).
Metropolitan Water Reclamation District of Greater Chicago (MWRDGC), 2008.
Environmental Monitoring and Research Division – 2007 Annual Report, Report No.
08-67.
Nevin, J.P., Connor, J.A., Newell, C.J., Gustafson, J.B., Lyons, K.A., 1997. Fate 5: A
Natural Attenuation Calibration Tool for Groundwater Fate and Transport Modeling,
NGWA Petroleum Hydrocarbons Conference, Houston, Texas, November 1997.
Pepper, I.L., Brusseau, M.L., Prevatt, F.J., Escovar, B.A., 2021. Incidence of PFAS in soil
following long-term application of class B biosolids. Sci. Total Environ. 793, 148449.
Roehler, K., Haluska, A.A., Susset, B., Liu, B., Grathwohl, P., 2021. Long-term behavior of
PFAS in contaminated agricultural soils in Germany. J. Contam. Hydrol. 241,
103812.
Sepulvado, J.G., Blaine, A.C., Hundal, L.S., Higgins, C.P., 2011. Occurrence and fate ofperfluorochemicals in soil following the land application of municipal biosolids.
Environ. Sci. Technol. 45 (19), 8106–8112.
Sharma, B.M., Bharat, G.K., Tayal, S., Larsen, T., Becanova, J., Karaskova, P.,
Whitehead, P.G., Futter, M.N., Butterfield, D., Nizzetto, L., 2016. Environ. Pollut.
208 (Part B), 704–713.
Silliman, S.E., Berkowitz, B., van Šimůnek, J., Genuchten, M.Th., 2002. Fluid flow andsolute migration within the capillary fringe. Ground Water 40 (1), 76–84.
Silva, J.A.K., Martin, W.A., Johnson, J.L., McCray, J.E., 2019. Evaluating air-water and
NAPL-water interfacial adsorption and retention of perfluorocarboxylic acids withthe vadose zone. J. Contam. Hydrol. 223, 103472 https://doi.org/10.1016/j.
jconhyd.2019.03.004.
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.
Silva, J.A.K., Šimůnek, J., McCray, J.E., 2022. Comparison of methods to estimate airwater interfacial areas for evaluating PFAS transport in the vadose zone. J. Contam.
Hydrol. 247, 103904 https://doi.org/10.1016/j.jconhyd.2022.103984.

14