Science of the Total Environment 935 (2024) 173252

Contents lists available at ScienceDirect

Science of the Total Environment
journal homepage: www.elsevier.com/locate/scitotenv

Numerical modeling of PFAS movement through the vadose zone: Influenceof plant water uptake and soil organic carbon distribution
Barbara Jennifer Biesek a, *, Adam Szymkiewicz a, Jirka Šimůnek b, Anna Gumuła-Kawęcka a,
Beata Jaworska-Szulc aabGdańsk University of Technology, Faculty of Civil and Environmental Engineering, ul. Gabriela Narutowicza 11/12, 80-233 Gdańsk, Poland
Department of Environmental Sciences, University of California Riverside, Riverside, CA, USA

H I G H L I G H T S

G R A P H I C A L A B S T R A C T• Use of organic carbon data from top soil
layers overestimates PFAS travel time.
• Root zone slows down PFAS travel time
due to transpiration.
• Root zone increases PFAS concentra­
tions in soil water.
• Influence of root zone on PFAS move­
ment is more important in fine soils.
• Root water uptake leads to increase of
PFAS accumulation on Air Water Inter­face Area.

A R T I C L E I N F O

A B S T R A C T

Editor: Damià BarcelóIn this study, we investigated the effects of soil organic carbon (SOC) distribution and water uptake by plant roots
on PFAS movement in the vadose zone with a deep groundwater table under temperate, humid climate condi­
tions. Two series of numerical simulations were performed with the HYDRUS computer code, representing theleaching of historical PFOS contamination and the infiltration of water contaminated with PFOA, respectively.
We considered soil profiles with three distributions of SOC (no SOC, realistic SOC distribution decreasing withdepth, and uniform SOC equal to the content measured in topsoil), three root distributions (bare soil, grassland,
and forest), and three soil textures (sand, sandy loam, and loam). The SOC distribution had a profound impact on
the velocity of PFOS movement. The apparent retardation factor for realistic SOC distribution was twice as large
as for the scenario with no SOC and more than three times smaller than for the scenario with uniformly high SOC
content. We also showed that the root distribution in soil profoundly impacts the simulations of PFAS migration
through soil. Including the root zone significantly slows down the movement of PFAS, primarily due to increased
evapotranspiration and reduced downward water flux. Another effect of water uptake by plant roots is an in­crease of PFAS concentrations in soil water (evapo-concentration). The evapo-concentration and the slowdown of
PFAS movement due to root water uptake are more significant in fine-textured soils than in sand.

Keywords:PFAS
HYDRUS
Root water uptake
Soil organic carbon
Soil contamination
Clean water and sanitation* Corresponding author.
E-mail address: bb2016@protonmail.com (B.J. Biesek).
https://doi.org/10.1016/j.scitotenv.2024.173252
Received 1 February 2024; Received in revised form 12 May 2024; Accepted 12 May 2024
Available online 18 May 20240048-9697/© 2024 Elsevier B.V. All rights are reserved, including those for text and data mining, AI training, and similar technologies.

B.J. Biesek et al.

Science of the Total Environment 935 (2024) 173252

sorption model was also shown by Zeng et al. (2021). Moreover, Brus­seau et al. (2021) and Zeng et al. (2021) found kinetic effects in AWI
sorption negligible for a wide range of conditions considered in theirstudies. Different approaches were also employed to calculate the AWI
area, which is a key variable governing PFAS sorption. AWI area can becalculated from the water retention curve, based on the thermodynamic
theory. This method was applied e.g. by Silva et al. (2020, 2022). In theabsence of detailed data, Wallis et al. (2022) approximated AWI area asa linear function of water saturation. Laboratory measurements in soils
suggest a nonlinear relationship between AWI area and water satura­tion, with measured AWI area for a specific saturation generally larger
than the one obtained from thermodynamic method (e.g. Brusseau andGuo, 2021). A possible explanation for this discrepancy is the roughnessof soil grains (Brusseau, 2023). Several studies used AWI – water satu­
ration relationship obtained directly from measurements (e.g. Guo et al.,2020) or multiplied the AWI area obtained from the thermodynamic
approach by an empirical scaling factor (e.g. Zeng et al., 2021). The
sorption of PFAS at AWI modifies the air-water surface tension andcapillary pressure in soil pores, leading to the so-called surfactantinduced flow. This effect was studied by Guo et al. (2020), Brusseau
et al. (2021), Zeng and Guo (2021) and Vahedian et al. (2024).
PFAS sorption on solid phase is a complex process. It occurs on soilorganic carbon (SOC), metal oxides and clay minerals and is affected byfactors such as pH, cation exchange capacity and soil micropore volume(Li et al., 2018; Nguyen et al., 2020). In modeling studies sorption on
solids is described using equilibrium models with linear isotherm (e.g.
Wallis et al., 2022) or Freundlich isotherm (e.g. Zeng et al., 2021; Silva
et al., 2022) or using kinetic models (e.g. Zeng et al., 2021). Zeng et al.
(2021) found kinetic effects in solid sorption to be of minor importancefor the field-relevant conditions examined in their study. Despite thecomplexity of PFAS sorption on solids, several authors used a simple
linear model, which assumes a direct proportionality between thepartition coefficient and the fraction of organic carbon in soils (e.g. Xiao
et al., 2015; Wallis et al., 2022; Smith et al., 2024). This model can beeasily parametrized and remains an attractive tool especially forscreening calculations (e.g. Smith et al., 2024). However, the SOC dis­tribution varies with depth. Typically, the maximum SOC concentration
occurs near the surface and becomes very small below 2–3 m depth(Jobbágy and Jackson, 2000). SOC is often estimated in samples takenfrom the shallow soil layer. For example, Xiao et al. (2015) used the SOC
fraction from shallow soil layers to estimate retardation factors for acontaminated site in Minnesota. Anderson (2021) noted that extrapo­
lating estimates of SOC fraction from the shallow horizons through awhole soil profile may lead to unrealistic results. Only few studies
accounted for the variability of solid sorption parameters with depth.
Wallis et al. (2022) assumed a variable SOC distribution was assumed,
based on soil sampling up to 7 m depth. The measured SOC profile was ingood agreement with the average profiles reported by Jobbágy andJackson (2000). Wallis et al. (2022) compared simulations withnonuniform SOC to simulations with no solid sorption and showed thatsorption on solids leads to strong retention of PFAS within the shallow
soil layer. Zeng and Guo (2023) used an extended linear sorption model,
where the partitioning coefficient depended on both SOC and the con­tent of clay and silt fraction (representing sorption on clay minerals).
They assumed SOC distribution decreasing exponentially with depth,
similar as in Jobbágy and Jackson (2000), which was constrained by themeasured SOC at multiple depths. Smith et al. (2024) assumed a ho­mogeneous vadose zone to apply analytical solutions, but the authors
used a thickness-weighted average of the soil organic carbon content,
which accounted for the fact that only a small fraction of SOC is present
below the depth of 0.6 m. However, the potential error in PFAS move­
ment prediction due to using SOC estimates from the top layer has notbeen quantified so far.
The relative importance of AWI and solid sorption for PFAS retention
in soil varies with the soil texture, heterogeneity, solid sorption prop­
erties and climatic conditions. For example, Guo et al. (2020) presented

1. Introduction
The term PFAS refers to a group of over 15,000 per- and poly­
fluoroalkyl substances containing carbon fluorine bonds (US EPA,
2024). Large amounts of energy are required to break a single C–F
bond, making it one of the strongest bonds in organic chemistry (Mur­phy et al., 2003). Possessing multiple C–F bonds is one of the charac­
teristics that give PFAS their stability. The widespread use of PFAS inmany industrial products and their resistance to environmental degra­
dation make them one of the most important groups of contaminants
worldwide. Although typically occurring in very low concentrations,they pose a considerable threat to human health. Increasing efforts are
directed at monitoring PFAS presence and investigating their behavior
in the environment. The soil compartment is particularly importantsince PFAS contamination typically originates at the ground surface, e.
g., from the application of aqueous film-forming foams (AFFFs) for fire
suppression, applications of biosolids, surface coating, de-icing compo­
nents, and exposure to contaminated industrial waste (Sharifan et al.,2021; Minnesota Pollution Control Agency, 2024). PFAS are sorbed on
soil particles and the air-water interface (AWI), potentially becoming along-term groundwater contamination source. PFAS fate in soils is also
affected by other factors, such as interactions with the microbiome andplant uptake (Sharifan et al., 2021).
Numerical models can be used to predict the fate of PFAS in soils andidentify the main factors that control it (e.g., Zeng et al., 2021; Silva
et al., 2022; Wallis et al., 2022). However, PFAS sorption on AWI de­pends on soil water content and cannot be represented by computercodes based on the standard transport equation (which typically in­cludes advection, dispersion, and different models of sorption on solid).
Until now, only a few numerical codes with the implementation of the
AWI sorption mechanism have been reported in the literature and are
publicly available. The well-known computer program HYDRUS(Šimůnek et al., 2008) based on the finite element method has been
modified to include sorption on AWI, first in 1D and 2D (Silva et al.,2020) and later in the full 3D version of HYDRUS 5 (Šimůnek et al.,2023). Another numerical model was developed by Guo et al. (2020) for
one dimension and later extended to three dimensions (Zeng and Guo,2021). In this computer program equations describing water flow and
PFAS transport are discretized using the cell-centered finite difference
method and fully implicit first-order discretization in time. The resulting
nonlinear equations are solved using Newton's iterative method. Wallis
et al. (2022) added a simplified representation of AWI sorption to the
LEACHM code (Hutson, 2003) describing 1D flow and transport in thesoil profile. Furthermore, simplified analytical solutions for PFAS
transport in the vadose zone were developed by Guo et al. (2022) andSmith et al. (2024). The above mentioned numerical and analyticalmodels were applied to investigate the role of various factors influencing
PFAS fate in vadose zone. Most applications focused on 1D verticaltransport in laboratory columns (Silva et al., 2020; Zeng et al., 2021), inmodel (synthetic) soil profiles (e.g. Silva et al., 2020; Guo et al., 2020;
Zeng et al., 2021) or in real soil profiles for which sufficient data wereavailable (Silva et al., 2022; Wallis et al., 2022). So far, only few studies
considered multidimensional flow (Zeng and Guo, 2021, 2023).
The AWI sorption of PFAS was described by equilibrium sorption
isotherms: linear (e.g. Wallis et al., 2022) or Langmuir (e.g. Silva et al.,2020; Vahedian et al., 2024), as well as by kinetic sorption models. Guo
et al. (2020) proposed a two-domain kinetic model for AWI sorption,
where the equilibrium part was described by Langmuir isotherm, while
the non-equilibrium part was described by a first-order kinetic model.
This approach was followed in other studies (e.g. Zeng et al., 2021).
Brusseau et al. (2021) concluded that linear equilibrium sorption models
are good approximations if PFAS concentrations in water are 3 to 10
times lower than their critical concentrations, specific to each PFAS
compound. The critical concentrations for two common PFAS com­pounds: PFOS and PFOA were reported as 1 mg/L and 10 mg/L,
respectively (Brusseau et al., 2021). The adequacy of linear equilibrium
2

B.J. Biesek et al.

Science of the Total Environment 935 (2024) 173252

simulations of hypothetical PFAS transport scenarios in two soil types
under typical humid and semi-arid climate. In their simulations, AWI
sorption was the main mechanism for PFAS retention in the vadose zone,
leading to retardation factors of the order from several hundred up to
more than one thousand. The retention of PFAS at AWI was more sig­
nificant in sand than in a soil with higher content of fine fractions. Silva
et al. (2022) simulated PFAS leaching due to biosolids applications at
two test sites in Illinois (US). In the considered scenarios, the effect of
AWI sorption was much smaller than the sorption on solids due torelatively large water contents (and thus small air-water interfacial
areas) in fine-textured soils. The same study also showed that AWI
sorption was more important in a homogeneous silty sand profile than ina heterogeneous profile with several layers of fine-textured soils un­derlain by medium sand. This was consistent with earlier findings ofSilva et al. (2020), who found that textural heterogeneities in 1D soil
profile lead to faster downward movement of PFAS. In the studies ofmultidimensional transport by Zeng and Guo (2021, 2023) local textural
heterogeneities were found to greatly impact PFAS movement, often
leading to preferential transport much faster than predicted for homo­
geneous soil. Zeng and Guo (2023) identified reduction of the available
AWI due to increasing water content in preferential flowpaths as themain mechanism accelerating movement of long-chain PFAS, such as
PFOS, in heterogeneous soil.
Zeng et al. (2021) investigated PFAS leaching in hypothetical 4 m
deep soil profiles in humid and semi-arid climates using realistic sets ofdaily weather data. They performed a comprehensive study of the model
sensitivity with respect to several parameters. It was shown that the rankof the most sensitive parameters depends on the relative importance ofsolid-phase sorption vs. AWI sorption, which in turn depends on the
PFAS chain length and applied concentrations. For two common com­pounds, PFOA and PFOS, the most sensitive parameters included dis­persivity, the parameter alpha in the van Genuchten retention function(which additionally influences the calculated AWI area), saturated watercontent, parameters describing AWI properties, and the exponent in the
Freundlich isotherm representing sorption on solids.
Until now, relatively little attention has been paid to the role of water
uptake by plant roots in the transport of PFAS. The presence of roots
modifies flow velocities and water contents, and consequently also AWI
areas in the soil profile. Wallis et al. (2022) seem to be the only study
that included the root zone in a PFAS transport model. They simulated
PFAS transport in a 7 m deep vadose zone profile in Australia, based ondata obtained from a site affected by historical AFFF application. Wallis
et al. (2022) considered seasonal growth of grass cover at the site, withvarying root depth (up to 30 cm) and including the effect of plant growth
on evapotranspiration. They showed a significant role of PFAS evapoconcentration in the upper part of the soil profile, i.e., an increase in
PFAS concentrations due to evapotranspiration of water from the soil
(their model did not include PFAS uptake by plants). Evapoconcentration was suggested as the main mechanism leading to pro­longed retention of PFAS in the top 1–2 m of soil.
The main aim of our study is to systematically analyze the effects ofroot water uptake and nonuniform carbon distribution on PFAS move­
ment in the vadose zone in a temperate humid climate. We performed
numerical simulations for two contamination scenarios: leaching of
historical contamination and continuous infiltration of contaminatedrainwater, with data partly based on the contaminated site in Minnesota
described by Xiao et al. (2015). For each scenario, we considereddifferent combinations of root distribution, soil organic carbon distri­bution, and soil texture.

(Version 5) computer program (Šimůnek et al., 2023, 2024). We
assumed that the vertical water flow in the vadose zone can be describedby the 1D Richards equation:[()]∂θ(h) ∂
∂h
=
K(h)
+1
− S
(1)
∂z
∂t
∂z
where θ is the volumetric water content (− ), t is time (T), z is the verticalcoordinate (oriented positively upwards) (L), K is the soil hydraulicconductivity (L T− 1), h is the water pressure head (negative in unsatu­
rated conditions) (L), and S is water uptake by plant roots (T− 1) (L, T,
and M denote the physical dimension of length, time, and mass
throughout the text). The soil hydraulic properties are described by thevan Genuchten - Mualem model (van Genuchten, 1980; Mualem, 1976):
θ(h) = θr + (θs − θr ) Se (h) = θr + (θs − θr )[1 + |αh|n ]K(h) = Ks Kr (Se (h) ) = Ks

− m√̅̅̅̅̅[)m ]2(
Se 1 − 1 − S1/m
e

(2)
(3)

where Se is the effective saturation (− ), θr is the residual water content(− ), θs is the water content at saturation (− ), α (L− 1), n (− ), m (− ) areparameters related to the pore size distribution (m = 1–1/n), and Ks is
the hydraulic conductivity at saturation (L T− 1). Simulations were per­formed for three soil textures: sand, sandy loam, and loam, using their
typical parameters taken from Carsel and Parrish (1988) listed inTable 1.
Water uptake by plant roots is described as:
S = a(h)b(z)Tp(4)

where a is the water stress response function for plant water uptake (− ),b is the normalized root density function (L − 1), and Tp is the rate ofpotential transpiration (L T− 1) (described later). The stress response
function was taken according to the well-known model of Feddes et al.
(1978). According to this model a = 0 for h > h1 and h < h4, which
represents conditions too wet or too dry for plants to extract water. The
optimum extraction rate (a = 1) occurs between h2 and h3 (h1 < h2 < h3< h4). Linear interpolation of a is applied between h1 and h2 and between
h3 and h4. In our study, we assumed h1 = − 10 cm, h2 = − 25 cm, h3 =− 800 cm, and h4 = − 8000 cm, while h3 varied according to the tran­
spiration rate, with h3 = − 200 cm for the rate of 0.5 cm day− 1 and h3 =− 800 cm for the rate of 0.1 cm day− 1. These are the default parameters
for pasture in HYDRUS 5, based on Wesseling et al. (1991).
The root distribution was taken after Jackson et al. (1996), who
proposed a formula describing the average distribution of roots withdepth in different biomes based on a compilation of existing literaturedata:
Y(z) = 1 − β|z|(5)where Y(z) is the cumulative root fraction (0 ≤ Y(z) ≤ 1) between thesoil surface and depth z (given as a positive number in cm) and β is aparameter depending on the type of biome. We performed simulations
for grass (temperate grassland, β = 0.943) and forest (temperate conif­
erous forest, β = 0.976). The corresponding cumulative root distribu­
tions are shown in Fig. 1. In both cases, the cumulative root distribution
was converted to the root density function b(z) from Eq. (4) by calcu­
lating its derivative and normalizing the result so that the integral of b(z)Table 1
Hydraulic parameters for three soil textures (from Carsel and Parrish 1988).

2. Materials and methods
2.1. Modeling water flow and root water uptake
All numerical simulations were performed with the HYDRUS

3Parameter

Sand

Sandy loamLoamθs (− )θr (− )
Ks (cm day− 1)
α (cm− 1)n (− )

0.43
0.045
712.8
0.145
2.68

0.41
0.065
106.1
0.075
1.89

0.43
0.078
24.96
0.036
1.56

B.J. Biesek et al.

Science of the Total Environment 935 (2024) 173252

zone can be described by the following equation (Silva et al., 2020; Guoet al., 2020):()∂(θ Cw ) ∂(ρb Cs ) ∂(Aaw Γ) ∂∂Cw∂(q Cw )
+
+
=
θD
−
(7)
∂z
∂t
∂t
∂t
∂z
∂z
where Cw is the concentration of PFAS in pore water (M L− 3), Cs is theconcentration of PFAS sorbed on solid phase (− ), ρb is the dry bulk
density of soil (M L− 3), Γ is the concentration of PFAS sorbed on AWI (M
L− 2), Aaw is the AWI interfacial area per unit soil volume (L− 1), D is thehydrodynamic dispersion coefficient (L2 T− 1), and q is the volumetricwater flux (Darcy velocity) (L T− 1). The hydrodynamic dispersion co­efficient is calculated as D = αL (q/θ) + τ Dm, where αL is the longitudinal
dispersivity [L], q/θ is the pore water velocity, τ is the tortuosity factor
for diffusion (− ) (τ = θ 7/3 / θs), and Dm is the molecular diffusion co­efficient for solutes in free water (L2 T− 1) (Table 1). The three terms on
the left-hand side represent the change of PFAS mass dissolved in water,
sorbed on the solid phase, and sorbed on the AWI, respectively.
Uptake of PFAS by plants was not considered in this study. We
considered only water uptake through plant roots to examine the effectsof water velocity and water content changes on PFAS movement in thesoil profile.
We assumed that the solid phase concentration is related to theconcentration in water by a linear isotherm:

Fig. 1. Distribution of plant roots with depth for grassland and forest (afterJackson et al., 1996).

Cs = Kd Cw = foc Koc Cw

over the soil profile is equal to 1. Note that Y(z) < 1 for all z, i.e. thenormalized root density is non-zero everywhere in the profile. However,the cumulative root fraction is equal to 0.9999 at depth of 158 cm forgrass and 380 cm for temperate coniferous forest. In our simulations we
set the relative density equal to zero below these depths. For compari­son, we also performed simulations for bare soil (without roots), as
described later.
In all simulations, the bottom boundary condition represented thegroundwater table, constant in time (h = 0). The profile depth was 2360
cm in the first series of simulations and 1372 cm in the second series ofsimulations. An atmospheric boundary condition with instantaneous
runoff (no ponding) was applied at the soil surface. Depending on thedaily weather data, the surface flux switched between infiltration (equal
to precipitation) and potential evaporation (Ep), defined as a fraction ofpotential evapotranspiration (PET). PET was assumed to be equal to thereference evapotranspiration calculated from the Penman-Monteith
formula (Allen et al., 1998). PET was split between Ep and TP (poten­
tial transpiration) according to de Beer's law, as implemented in HYD­
RUS (Šimůnek et al., 2023):
Ep = PET exp( − k • LAI) Tp = PET (1 − exp( − k • LAI) )(8)

where Kd is the linear distribution coefficient (L3 M− 1), foc is the fractionof organic carbon in soil (− ), and Koc is the partitioning coefficient for a
substance in pure organic carbon (L3 M− 1). Koc values are typicallydetermined through laboratory measurements or referenced fromavailable literature. The values used in our study are reported in Table 2.
The dry bulk density of soil was assumed to be 1.59 g cm− 3 in all sim­
ulations, following the study of Xiao et al. (2015).
The fraction of organic carbon in soil decreases with depth. Average
distributions of soil carbon in the top 3 m of soil were reported by
Jobbágy and Jackson (2000) for different biomes. The carbon distribu­
tions were provided as a fraction of the total mass of soil carbon in aninterval of known soil depths. The intervals allowed for the estimation of
foc distribution in known depths if the value of foc in the shallow soil
layer (0–20 cm) was known. In order to construct a soil profile with arealistic distribution of carbon, we assumed that foc in the top soil layer
up to a 20 cm depth was equal to 1.25 %, as reported by Xiao et al.
(2015). We used the OC distribution from Jobbágy and Jackson (2000)
to calculate foc at larger depths.
The distribution of foc and Kd to the 3 m depth is shown in Table 3,
which corresponds to data distribution for temperate grasslands from
Jobbágy and Jackson's (2000) study. The SOC distributions reported fortemperate forests and bare soils (desert) showed minor differences. For
the sake of simplicity and focus on the key priorities of this research, weused the same SOC with all three types of plant cover. Due to a lack ofavailable data for foc and Kd below the 300 cm depth, we assumed they
were constant and the same as in the 200–300 cm interval. For com­
parison, we also ran a simulation, where we assumed a complete lack of
SOC (and sorption only occurred to the AWI) and a constant Kd value(corresponding to the value used throughout the Xiao et al. (2015)
study), corresponding to the one from the topmost soil layer (0–20 cm).
Brusseau et al. (2021) demonstrated that depending on measured
surface tension and for the PFAS concentrations 3 to 10 times lower than(6)

where LAI is the leaf area index (− ), and k is the radiation extinction
coefficient, taken as 0.463. In the bare soil scenarios, LAI = 0 and Ep =PET; thus, the whole potential evapotranspiration flux was assigned as
the boundary condition at the surface. For both grass and forest, we
assumed LAI = 2.88, corresponding to the reference evapotranspiration
conditions. In general, LAI depends on plant species and changes during
the year. LAI influences the value of PET calculated from PenmanMonteith, as well as its distribution between Ep and Tp via de Beer's
law. However, our earlier study (Gumuła-Kawęcka et al., 2022) obtained
satisfactory modeling results using LAI values from 2.1 to 2.9 forgrassland and a pine forest in a region with a similar climate. Whileusing the same constant value of LAI for both types of plant cover is a
simplification of reality, it allows the effects of root distribution on waterflow to be separated from the impact of varying potential transpiration.

Table 2
Parameters describing PFAS sorption and diffusion (from Zeng et al., 2021).

2.2. Modeling PFAS transport
In this study, we considered two common compounds from the PFAS
group: PFOA and PFOS. The movement of each compound in the vadose
4

Parameter

PFOA

PFOS

Koc (cm3 g− 1)
Kaw (cm)
Dm (cm2 d− 1)105
0.00580.467

617
0.04270.467

B.J. Biesek et al.

Science of the Total Environment 935 (2024) 173252

simulations. The initial condition for Eq. (1) was obtained by running awarm-up simulation for 100 years, starting from the hydrostatic pres­sure head distribution above the water table. The final distribution of thewater pressure head was used as the initial condition in the actual
simulation of the PFOS transport. The initial distribution of PFOS
approximated the one shown in Fig. 3 of the paper by Xiao et al. (2015),
labeled Site 27. The total PFOS concentration in the soil increased lin­
early from 0 at the surface to the peak value of 160 ng per gram dryweight (equal to 0.254 μg per cm3 of soil) at a depth of 68 cm and thendecreased linearly to 0 at a depth of 90 cm. In all simulations, the lon­
gitudinal dispersivity αL was equal to 24 cm. This value was used in the
study by Zeng et al. (2021), although for a shorter soil profile. Zeng et al.
(2021) showed that the choice of αL significantly influences PFAS
movement in the subsurface. However, in this simulation series, we
focus on the movement of peak concentration, which is affected mostly
by the advection velocity and the retardation factor.
We considered three types of root distributions: bare soil, grass, andforest, and three types of SOC distributions: SOC-0 (no SOC in the soil
profile), SOC-V (variable distribution of SOC, decreasing with depth, as
described above), and SOC-H (a uniformly high value of SOC, corre­
sponding to the one from the top soil layer, i.e., 1.25 %, as used by Xiao
et al., 2015). Altogether, nine simulations were performed to represent
PFOS movement for all the possible combinations of the above factors.

Table 3
Nonuniform distribution of organic carbon in a soil profile.
Depth [cm]% of the total SOC storage*

foc [%]

0–20
20–40
40–60
60–80
80–100100–200
200–300

41
23
15
12
9
7
5

1.250
0.701
0.457
0.366
0.274
0.219
0.167

*

Kd (cm3 g− 1)
PFOA

PFOS

1.309
0.734
0.479
0.383
0.287
0.229
0.175

7.708
4.324
2.820
2.256
1.692
1.350
1.028

Within the top 3 m of soil.

their critical reference concentrations (as measured by Brusseau and
Van Glubt, 2019), the adsorption of contaminants to the air-waterinterfacial area (AWI) follows a linear process. According to this linear
equilibrium model, the PFAS concentration at AWI, Γ, can be deter­
mined as:
Γ = Kaw Cw(9)where Kaw is the coefficient for AWI adsorption in [cm3 cm− 2]. The
coefficients for PFOA and PFOS used in our study are given in Table 2.
The AWI area is calculated as:∫ρ g θsAaw = wh(Θ)dΘ
(10)

σ02.4. Simulation setup for infiltration of water contaminated with PFASθw

In the second simulation series, we investigated the effect of root
distribution and soil texture on the movement of PFOA. We considered ahomogeneous soil profile with a groundwater table at a depth of 1372
cm. The profile was discretized using 937 nodes, with node spacing of 1
cm in the upper 500 cm of the profile and 2 cm in the lower part.
Simulations were performed for 9 cases, i.e., three soil textures: sand
(Sa), sandy loam (SaLo), and loam (Lo), and three types of root distri­butions: bare soil, grass, and forest, as described previously. We used thesame realistic SOC distribution in all simulations, with the OC content
decreasing with depth (SOC-V), as described earlier. We used daily
weather data from January 1, 1979, to July 31, 2014 (12,987 days),taken from the same source as in the first simulation series. In the case of
the forest on the loam soil, we had to extend the simulation time byrepeating weather data from the earlier period due to the very slow
movement of PFOA. The infiltrating water has a constant PFOA con­centration of 1 μg cm− 3.
At the bottom of the profile, the boundary conditions were h = 0 forwater flow (a groundwater table, constant in time) and a zero concen­
tration gradient (for PFOA transport). A water flow simulation, startingfrom the hydrostatic pressure head distribution, was performed using
the daily weather data for 35 years to obtain realistic initial conditions
for water flow. The final pressure head profile was taken as the starting
condition for the actual simulation. The initial PFOA concentration was
set to 0 in the whole profile. The longitudinal dispersivity αL was equal
to 63 cm in all simulations (estimated using the formula of SchulzeMakuch, 2005). We compared the results in terms of PFOA concentra­
tions in water, water fluxes, and the cumulative PFOA mass that passed
through the depth of 2 m and 6 m in the soil profile.

where ρw is the water density (M L− 3), g is the gravitational acceleration(L T− 2), σ0 is the air-water surface tension (M T− 2), and Θ is an inte­
gration variable (− ). This approach is implemented in HYDRUS 5,
following Bradford and Leij (1997) and Bradford et al. (2015). Aaw in­creases with decreasing water content.
The retardation coefficient R (− ) is commonly used to quantify how
much the movement of a sorbing contaminant is slowed down compared
to a conservative (non-sorbing) contaminant. Assuming linear sorption
on the solid phase and AWI, the retardation coefficient can be defined as:ρ Kd Aaw Kaw
R = 1+ b +
θ
θ

(11)

where R can be meaningfully defined when θ, Aaw, and the specificdischarge q are approximately constant in space and time, i.e., for quasisteady, uniform flow.
2.3. Simulation setup for leaching of historical PFAS contamination
The first series of simulations was partly based on the data from acontaminated site in Minnesota (USA), described by Xiao et al. (2015).
These simulations aimed to show the effects of different assumptions
about SOC and root distribution in soil. Only the PFOS movement wasconsidered. The soil profile was homogeneous and consisted of sandyloam. The profile depth was 2360 cm, based on the depth of the watertable in the location corresponding to Well #1, which was close to Site
27 and 28, where Xiao et al. (2015) obtained groundwater samples fortheir study. The profile was uniformly discretized with 937 nodes. The
weather data source was the Climate Forecast System Reanalysis (CFSR),
accessible as part of the Soil & Water Assessment Tool (SWAT) at htt
ps://swat.tamu.edu/data/cfsr. We considered daily weather data from1979 to 2013 and chose a representative year with precipitation andpotential evapotranspiration (PET) close to the average annual values
calculated for the entire 35-year period. Out of the 35-year period, thedata from 2002 was used since the average potential evapotranspiration(1024 mm/year) and the annual precipitation (903 mm/year) was
closest to the average of the 25 years (1022 mm/year and 851 mm/year
respectively). The data for the selected year were then repeated 100
times to create a 100-year series of typical weather data used in the

3. Results and discussion
3.1. Leaching of historical PFAS contamination
Fig. 2 shows the PFOS concentrations in the soil profile after 100
years for all nine scenarios. The position and magnitude of the peakconcentration are strongly influenced by the assumed root distribution
and SOC distribution. The fastest movement and smallest peak attenu­
ation occur in bare soil, followed by grass and forest. The presence of
SOC significantly slows down the movement of PFOS and decreases its
peak concentration. Predictably, this retardation effect is much stronger
5B.J. Biesek et al.

Science of the Total Environment 935 (2024) 173252

Fig. 2. Profiles of PFOS concentrations in water after 100 years of simulation time in bare soil (left), grassland (middle), and forest (right) (the historicalcontamination leaching scenario).

topsoil (the same as in our study) and Koc = 479 cm3 g− 1 (less than in ourstudy). Moreover, they neglected AWI sorption in the calculation of R.
Our estimate of R with a realistic SOC distribution varies between 29 and
32 (Table 1), about two times less than in Xiao et al. (2015), despite thepresence of AWI sorption. Without AWI sorption, the R factor for the
SOC-V scenario would be only about 11 to 13 (the difference between R
for SOC-V and SOC-0 in Table 4 plus one). Based on our simulation re­sults, for SOC-H without AWI sorption, one would obtain R around 78 to
88, somewhat larger than 64 reported by Xiao et al. (2015), due to thelarger Koc used in our study.
The retardation factors calculated from the travel distance of thepeak concentration in the numerical solution are in good agreementwith the estimates from Eq. (11) for those scenarios where foc and Kd are
uniform, i.e., for SOC-0 and SOC-H. The evolution of the AWI area
within a single year for different scenarios is shown in Fig. 3. Note thatthe same pattern repeats each year because of the repeating set of yearly
weather data. Significant temporal changes in the AWI area are visible
only in the scenarios with the root zone in the uppermost part of the
profile, where θ decreases and Aaw increases due to increased evapo­transpiration in summer. Taking approximately Aaw = 70 cm− 1, which isrepresentative of the deeper vadose zone, we obtain from Eq. (11) R =20.9 for θ = 0.15 and Kd = 0 (SOC-0) and R = 102.5 for Kd = 7.7 cm3 g− 1
(SOC-H). These estimates are reasonably close to the values in Table 1
because the parameters involved in Eq. (1) do not change significantly inspace and time. However, estimating Kd that should be used in Eq. (1) for
the SOC-V scenario is challenging, as Kd varies with depth. The average
Kd can be back-calculated from Eq. (1), assuming that (without AWI) R= 1 + ρb Kd/θ should be approximately equal to 12 (see above). The
resulting Kd is 1.04 cm3 g− 1, i.e., very close to the minimum Kd value
applied in the soil profile below the depth of 200 cm. Thus, the long-in the SOC-H scenarios than in SOC-V. It can be noted that the results for
the SOC-V scenarios are quite different from both SOC-H and SOC0 scenarios.
The distance traveled by the peak concentration in each scenario
from the initial position of the peak at a depth of 68 cm is reported inTable 4. For comparison, we also estimate the distance traveled in thesame time frame by a conservative (non-sorbing) contaminant. It is
calculated as L = (q/θ) t, where q is the average unit discharge throughthe soil profile, θ is the average water content, and t = 100 years. The
average unit discharge was assumed to equal the average flux reaching
the water table at the bottom of the soil profile (groundwater recharge).
It equals 42.2 cm yr− 1 for bare soil, 24.3 cm yr− 1 for grassland, and 19.0
cm yr− 1 for forest. Fig. 3 shows the water content changes theta in both
space and time. The flow appears quasi-steady and θ does not vary much
between the scenarios. We assumed θ = 0.15 to calculate the distances
for the conservative contaminant shown in Table 4. This value is in good
agreement with Fig. 3. Coincidentally, it is also equal to the value usedby Xiao et al. (2015) to estimate the retardation factor in their study,
which provided the basis for our numerical simulations.
In Table 4, we also present retardation factors estimated as R = L0/L,
where L0 is the distance traveled by the conservative contaminant, and L
is the distance traveled by PFOS for the same root distribution. R differs
significantly depending on the SOC distribution. Even for the SOC0 scenario, we obtained significant retardation factors around 18–19.
This is due to the sorption on AWI, which occurs regardless of theassumed SOC distribution. The retardation factor for the SOC-V scenario
is more than three times smaller than for SOC-H, which emphasizes theneed for realistic data on SOC to predict PFAS fate in the vadose zone.
For instance, Xiao et al. (2015) estimated the retardation factor for PFOS
as R = 64, assuming ρb = 1.59 g/cm3, θ = 0.15, foc = 1.25 % in the

Table 4
Travel distances of the peak concentration in 100 years and retardation factors calculated for different scenarios of historical contamination leaching.
Scenario
Conservative
SOC-0
SOC-V
SOC-H

Bare soil

Grass

ForestDistance (cm)R (− )

Distance (cm)R (− )

Distance (cm)R (− )

28,122
1523
930
275

1.00
18.46
30.24
102.26

16,196
865
509
151

1.00
18.72
31.82
107.26

12,685
676
426
131

1.00
18.76
29.78
96.83

6

B.J. Biesek et al.

Science of the Total Environment 935 (2024) 173252

Fig. 3. Evolution of the volumetric water content (top) and AWI area (bottom) in bare soil (left), grassland (middle), and forest (right) during a single year in thesimulation of historical PFOS contamination leaching.

term fate of PFAS in the deep vadose zone is not significantly influenced
by the increased organic carbon content of the soil close to the surface.
We should note that in real settings, the carbon content in the deep
vadose zone is probably even lower than 0.167 %, as used here. Again,
this shows that realistic estimates of SOC in deeper soil horizons are ofparamount importance for PFAS modeling.
The assumptions regarding root distribution also significantly in­fluence the predictions of PFOS movement. Fig. 3 shows that the pres­ence of roots clearly increases the AWI area in the root zone because thewater content is reduced due to water uptake by plants. Thus, one could
expect enhanced sorption of PFOS on AWI in the scenarios with plant
cover. However, as shown in Table 4, all retardation factors, including
those related to AWI sorption only (the SOC-0 scenario), are quite
similar for each type of root distribution. It seems that, similarly to theincreased SOC, the larger Aaw due to transpiration in the upper soil
layers does not significantly affect the long-term movement of PFAS. The
main mechanism by which the roots influence PFAS movement seems tobe the reduction in the downward water flux due to increased evapo­transpiration. Using data from Table 1, we can calculate an apparent
retardation factor due to the presence of root zone as Rr = LB / L, where
LB is the distance traveled by the peak concentration in the bare soil
profile, and L is the corresponding distance in a profile with vegetation.
Rr was found to be quite similar for each SOC distribution, ranging from1.76 to 1.83 for grass and from 2.10 to 2.25 for forest. These numbers
closely match the corresponding ratios of the groundwater recharge
rates, which are 1.74 (bare soil/grass) and 2.22 (bare soil/forest). A
similar effect of slowing the downward movement due to the root zone
was also shown for conservative contaminants (e.g., Szymkiewicz et al.,2018, 2019).

show oscillations, with magnitudes depending on the soil type and root
distribution. We note that these changes are related to the dynamic
nature of the AWI area, which dramatically changes with changes inwater contents under transient water flow conditions, dynamicallyaltering the AWI sorption capacity. This is contrary to sorption to thesolid phase, which is constant (proportional to the soil bulk density) andindependent of transient flow conditions. The oscillations in the sand are
relatively small, and quasi-steady conditions are reached at both 2 m
and 6 m depths. The average concentration in the quasi-steady phase islarger than in the infiltrating water: around 1.4 mg/cm3 for bare soil,
2.2 mg/cm3 for grass, and 2.7 mg/cm3 for forest. The PFOA concen­
tration in water increases due to evapotranspiration, which takes out
some of the water and reduces the downward flux without decreasing
the mass of dissolved PFAS. This effect is more pronounced for grass andforest than for bare soil since the presence of the root zone increases
actual evapotranspiration. Moreover, the increase in concentration ismore significant fine textured soils. In Fig. 4 the maximum PFOA con­centrations in pore water are higher in sandy loam (around 5 μg cm− 3)
and loam (around 7 μg/cm− 3) than in sand (around 3 μg/cm− 3). The
increase of PFAS concentration due to evapotranspiration can be called
evapo-concentration, following Wallis et al. (2022). However, our sim­
ulations were performed for a humid climate with general excess ofinfiltration over evapotranspiration and the upward water fluxes in soil
were negligible, as shown in Figs. 4 and 5. In contrast, the study of Wallis
et al. (2022) focused on a more arid climate, where the upward move­
ment of water with dissolved PFAS significantly contributed to increased
concentrations in top soil.
The oscillations in PFOA concentrations seen in Figs, 4 and 5 are
related to the seasonal changes in water contents, leading to changes in
AWI areas. In sandy loam and loam, the changes in water contents and
AWI are much larger than in sand (Fig. 6), causing larger oscillation
amplitudes in the PFOA concentrations. The maximum PFOA concen­
trations in loam reach up to 7 mg/cm3 for the forest. Also, the steadystate is not reached in loam even at the end of simulation time because

3.2. Infiltration of water contaminated with PFAS
The results of simulations are shown in Figs. 4 and 5 for the depths 2
and 6 m, respectively. The breakthrough curves for PFOA concentrations
7B.J. Biesek et al.

Science of the Total Environment 935 (2024) 173252

Fig. 4. PFOA concentrations in water (top), water fluxes (middle), and cumulative PFOS mass fluxes (bottom) at a depth of 2 m in sandy (left), sandy loam (middle)
and loam (right) soil profile in the simulations representing infiltration of contaminated water. Results for different root distributions are distinguished by colors(blue – bare soil, green – grass, red – forest).

the flow and transport processes are much slower than in sand or sandyloam. The decrease in water content during dry months lead to coun­teracting effects on PFAS concentrations. On the one hand, a decrease inwater content without a loss of PFOA and PFOS should lead to higher
concentration in the bulk soil. On the other hand, an increase in AWI
leads to lower concentration in the water phase.
The variability of weather conditions between subsequent years also
influences the results. For example, in Fig. 4, in the loam-forest scenario,
a significant decrease in concentrations is observed approximately be­
tween day 3000 and 4500. Fig. 4 shows that this is a dry period withalmost no flow in the soil profile. On the other hand, we can see in Fig. 6that in the corresponding time period, the water content is reduced, and
the AWI area increased in the part of the loam profile around 2 m deep,
leading to enhanced sorption of PFOA on AWI. Note that a different
color scale is used for the AWI area in Fig. 6 for the loam and sandprofiles. The maximum Aaw in sand during the whole simulation is67.11 cm− 1, while in loam, it is 2059 cm− 1.
Since concentrations and water fluxes vary significantly with time,an analysis of concentration evolution is insufficient to quantify theslowdown of contaminant movement due to the presence of roots. In thiscontext, the plots of the cumulative PFOA mass passing through a spe­cific depth are more helpful (shown in the bottom rows of Figs. 4 and 5).
While PFOA transport in profiles with root zones is always slower thanin the bare soil profile, the difference depends strongly on the soil
texture and depth and varies with time. For example, choosing

arbitrarily the characteristic time corresponding to the cumulative massequal to 500 μg, we obtain the ratios of times between forest and baresoils given in Table 5 (for loam, we had to extend the simulation period;
the results are not shown in Fig. 5). The ratios are clearly larger for finer
textures (sandy loam and loam) than for sand, and they increase withdepth for all soil textures.
3.3. Limitations of the analysis
The presented results were obtained using several simplifying as­
sumptions, which must be clearly understood. We focused on two
common PFAS compounds: PFOS and PFOA. However, although fielddata show that these compounds are the most frequently tested PFAS,
other species are present in contaminated zones (e.g., Johnson, 2022;
Batayi et al., 2021; Silva et al., 2022). Since the PFAS group expandedfrom 9000 to 15,000 compounds (US EPA, 2021) just during thisresearch, many more compounds should be considered when deter­mining the real soil and water concentrations of these pollutants.
Furthermore, we assumed a simplified description of PFAS sorption on
AWI and SOC using linear models. AWI area was calculated using thethermodynamic approach, which may underestimate the extent ofsorption, especially in dry conditions (Brusseau, 2023). Furthermore,
solid sorption was related only to SOC, while it may also occur on clay
minerals and metal oxides (e.g. Zeng and Guo, 2023). We are also awarethat PFAS could be absorbed by plants and adsorbed to plant roots.
8B.J. Biesek et al.

Science of the Total Environment 935 (2024) 173252

Fig. 5. PFOA concentrations in water (top), water fluxes (middle), and cumulative PFOS mass fluxes (bottom) at a depth of 6 m in sandy (left), sandy loam (middle)
and loam (right) soil profile in the simulations representing infiltration of contaminated water. Results for different root distributions are distinguished by colors(blue – bare soil, green – grass, red – forest).

However, only limited data on this factor are currently available andvary depending on the type of plant species and PFAS, as well as infil­
trating concentrations (Held and Reinhard, 2020; Felizeter et al., 2021;
Xu et al., 2022; Blaine et al., 2014). Therefore, we could not incorporate
this factor into our study due to insufficient data.
For simulation models using data from temperate grasslands andtemperate coniferous forests, the average constant for a variety of spe­
cies present in these biomes was used (Jackson et al., 1996), rather thanaccounting for the diversity of species present, which would be a more
realistic representation. Also, we did not consider seasonal or long-termroot growth or the LAI changes affecting evapotranspiration.
We considered only homogeneous soil profiles, while it has beenshown that soil heterogeneity influences PFAS movement. For instance,
Silva et al. (2022) found that the impact of AWI adsorption was moresignificant for PFAS transport within homogeneous soils than texturally
heterogeneous soils considered in their study. Zeng and Guo (2021,2023) used a 3D model to show the importance of preferential flow
paths and capillary barriers on PFAS spreading in the subsurface. They
concluded that soil heterogeneities are among the primary sources of
uncertainty for predicting PFAS leaching and retention in the vadosezone. The interplay between water uptake by roots and soil

heterogeneities was outside the scope of this study and remains to beinvestigated.
In both series of simulations, the groundwater table was relatively
deep and constant in time. Groundwater table fluctuations can be ex­pected to affect PFAS fate by creating or destroying AWI during thetransition from unsaturated to saturated conditions and vice versa. The
weather conditions used to run the first series of simulations were ob­tained by selecting a yearly data set from the 35-year period to depict the
average conditions. The limitation of using this data is that it neglects
climate change and fluctuating weather patterns such as wet or dry
periods. The second series of simulations showed the influence of yearto-year variations in weather conditions.
4. Conclusions
In this study, we investigated the effect of soil organic carbon dis­tribution and water uptake by plant roots on PFAS movement in thevadose zone with a deep groundwater table under temperate, humid
climate conditions. We found that both SOC distribution and root zone
significantly affect PFAS transport.
The apparent retardation factor for realistic SOC distribution was
9B.J. Biesek et al.

Science of the Total Environment 935 (2024) 173252

Fig. 6. AWI areas (top) and water contents (bottom) as a function of depth and time during the simulations representing infiltration of contaminated water in the
sandy (left) and loamy (right) forest profiles.

considering various soil textures and SOC distributions and focusing on adifferent climate.
Despite several simplifying assumptions, our study highlights theimportance of using realistic data about the distribution of soil organiccarbon and plant roots for accurate predictions of PFAS fate in thesubsurface. Further research is needed to investigate the role of addi­
tional factors such as soil heterogeneity, groundwater table fluctuations,and climate change.

Table 5
Ratios between the characteristic times for forest and bare soils. The charac­
teristic time corresponds to when the cumulative leaching mass at a particulardepth of 500 μg is reached, as shown in Figs. 4 and 5.
Soil texture

Depth 2 m

Depth 6 m

Sand
Sandy loam
Loam

1.13
1.32
1.34

1.28
1.58
1.54

CRediT authorship contribution statement

twice as large as for the scenario with no SOC (only AWI sorption) andmore than three times smaller than for the scenario with a uniformly
high SOC content, equal to the content in topsoil. Thus, realistic SOC
distribution must be considered when predicting PFAS fate in soils andgroundwater vulnerability to contamination. Using SOC estimates basedon sampling shallow soil layers overestimates PFAS travel times
significantly.
We also showed that the root distribution in soil profoundly impacts
the simulations of PFAS migration through soil. The presence of the rootzone slows the movement of PFAS significantly, primarily due toincreased evapotranspiration and reduced downward water flux.
Another effect of water uptake by plant roots is an increase in PFAS
concentrations in soil water, i.e. evapo-concentration, shown earlier by
Wallis et al. (2022) for a more arid climate. In our simulations, signifi­
cant evapo-concentration was observed in a temperate humid climate.
PFAS concentrations changed seasonally with varying soil water con­tents. They were also influenced by year-to-year variability in precipi­
tation. The evapo-concentration and the slowdown of PFAS movement
due to root water uptake are more significant in fine-textured soils thanin sand. Our work extended the analysis of Wallis et al. (2022) by

Barbara Jennifer Biesek: Writing – original draft, Visualization,
Investigation, Formal analysis, Data curation, Conceptualization. Adam
Szymkiewicz: Writing – review & editing, Supervision, Methodology,
Investigation, Conceptualization. Jirka Šimůnek: Writing – review &editing, Software. Anna Gumuła-Kawęcka: Investigation. Beata
Jaworska-Szulc: Investigation.
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
Data will be made available on request.

10

B.J. Biesek et al.

Science of the Total Environment 935 (2024) 173252

References

Murphy, C., Schaffrath, C., O’Hagan, D., 2003. Fluorinated natural products: thebiosynthesis of fluoroacetate and 4-fluorothreonine in Streptomyces cattleya.
Chemosphere 52, 455–461. https://doi.org/10.1016/S0045-6535(03)00191-7.
Nguyen, T.M.H., Bräunig, J., Thompson, K., Thompson, J., Kabiri, S., Navarro, D.A.,
Mueller, J.F., 2020. Influences of chemical properties, soil properties, and solution
pH on soil–water partitioning coefficients of per-and polyfluoroalkyl substances(PFASs). Environ. Sci. Technol. 54 (24), 15883–15892.
Schulze-Makuch, D., 2005. Longitudinal dispersivity data and implications for scaling
behavior. Groundwater 43 (3), 443–456.
Sharifan, H., Bagheri, M., Wang, D., Burken, J.G., Higgins, C.P., Liang, Y., Liu, J.,
Schaefer, C., Blotevogel, J., 2021. Fate and transport of per-and polyfluoroalkyl
substances (PFASs) in the vadose zone. Sci. Total Environ. 771, 145427.
Silva, J.A., Guelfo, J.L., Šimůnek, J., McCray, J.E., 2022. Simulated leaching of PFAS
from land-applied municipal biosolids at agricultural sites. J. Contam. Hydrol. 251,
104089.
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 (10), 2758.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2008. Development and applications of the
HYDRUS and STANMOD software packages and related codes. Vadose Zone J. 7 (2),587–600.
Šimůnek, J., van Genuchten, M.T., Šejna, M., 2023. The HYDRUS Software Package forSimulating the One-, Two-, and Three-Dimensional Movement of Water, Heat, andMultiple Solutes in Variably-Saturated Porous Media. Technical Manual II. HYDRUS
2D/3D. ver. 5.xx.
Šimůnek, J., Brunetti, G., Jacques, D., van Genuchten, M.Th., Šejna, M., 2024. Recent
developments and applications of the HYDRUS computer software packages since2016, a special issue tribute to Rien van Genuchten, recipient of the 2023 wolf prize
for agriculture. Vadose Zone J. https://doi.org/10.1002/vzj2.20310 accepted
January 2, 2024.
Smith, J., Brusseau, M.L., Guo, B., 2024. An integrated analytical modeling frameworkfor determining site-specific soil screening levels for PFAS. Water Res. 252, 121236.
Szymkiewicz, A., Gumuła-Kawęcka, A., Potrykus, D., Jaworska-Szulc, B., PruszkowskaCaceres, M., Gorczewska-Langner, W., 2018. Estimation of conservative contaminant
travel time through vadose zone based on transient and steady flow approaches.
Water 10 (10), 1417.
Szymkiewicz, A., Savard, J., Jaworska-Szulc, B., 2019. Numerical analysis of recharge
rates and contaminant travel time in layered unsaturated soils. Water 11 (3), 545.
United States Environmental Protection Agency, 2021. EPA 822-F-21-007. TechnicalFact Sheet: Human Health Toxicity Assessment for GenX Chemicals. Table 2.
Comparison of Chronic Toxicity Values for the Four PFAS with Final EPA Toxicity
Assessments. Retrieved August 19, 2022 from. https://www.epa.gov/system/files/
documents/2021-10/genx-final-tox-assess-tech-factsheet-2021.pdf.
United States Environmental Protection Agency, 2024. Comptox Chemicals Dashboard:
Master List of PFAS Substances (Version 2). Retrieved Jan 13, 2024 from. https://comptox.epa.gov/dashboard/chemical_lists/pfasmaster.
Vahedian, F., Silva, J.A.K., Šimůnek, J., McCray, J.E., 2024. The influence of tensiondriven flow on the transport of AFFF in unsaturated media. ACS ES&T Water 11.
https://doi.org/10.1021/acsestwater.3c00611.
Wallis, I., Hutson, J., Davis, G., Kookana, R., Rayner, J., Prommer, H., 2022. Model-based
identification of vadose zone controls on PFAS mobility under semi-arid climate
conditions. Water Res. 225, 119096.
Wesseling, J.G., Elbers, J.A., Kabat, P., van den Broek, B.J., 1991. SWATRE: Instructions
for Input, Internal Note. Winand Staring Centre, Wageningen, the Netherlands.
Xiao, F., Simcik, M.F., Halbach, T.R., Gulliver, J.S., 2015. Perfluorooctane sulfonate(PFOS) and perfluorooctanoate (PFOA) in soils and groundwater of a US
metropolitan area: migration and implications for human exposure. Water Res. 72,
64–74.
Xu, B., Qiu, W., Du, J., Wan, Z., Zhou, J.L., Chen, H., Liu, R., Magnuson, J.T., Zheng, C.,
2022. Translocation, bioaccumulation, and distribution of perfluoroalkyl andpolyfluoroalkyl substances (PFASs) in plants. iScience Mar 11;25 (4), 104061.
https://doi.org/10.1016/j.isci.2022.104061. PMID: 35345465; PMCID:
PMC8957016.
Zeng, J., Guo, B., 2021. Multidimensional simulation of PFAS transport and leaching in
the vadose zone: impact of surfactant-induced flow and subsurface heterogeneities.
Adv. Water Resour. 155, 104015.
Zeng, J., Guo, B., 2023. Reduced accessible air–water interfacial area accelerates PFAS
leaching in heterogeneous vadose zones. Geophys. Res. Lett. 50 (8),e2022GL102655.
Zeng, J., Brusseau, M.L., Guo, B., 2021. Model validation and analyses of parameter
sensitivity and uncertainty for modeling long-term retention and leaching of PFAS in
the vadose zone. J. Hydrol. 603, 127172.

Allen, R.G., Pereira, L.S., Raes, D., Smith, M., 1998. Crop evapotranspiration-guidelines
for computing crop water requirements-FAO irrigation and drainage paper 56. Fao,
Rome 300 (9), D05109.
Anderson, R.H., 2021. The case for direct measures of soil-to-groundwater contaminantmass discharge at AFFF-impacted sites. Environ. Sci. Technol. 55 (10), 6580–6583.
Batayi, B., Okonkwo, O., Daso, A.P., 2021. Poly- and perfluorinated substances inenvironmental water from the Hartbeespoort and Roodeplaat Dams, South Africa.
Water SA 47. https://doi.org/10.17159/wsa/2021.v47.i1.9445.
Blaine, A., Rich, C., Sedlacko, E., Hyland, K., Stushnoff, C., Dickenson, E., Higgins, C.,
2014. Perfluoroalkyl acid uptake in lettuce (Lactuca sativa) and strawberry (Fragaria
ananassa) irrigated with reclaimed water. Environ. Sci. Technol. 48 https://doi.org/
10.1021/es504150h.
Bradford, S.A., Leij, F.J., 1997. Estimating interfacial areas for multi-fluid soil systems.
J. Contam. Hydrol. 27, 83–106.
Bradford, S.A., Wang, Y., Torkzaban, S., Šimůnek, J., 2015. Modeling the release of E.
Coli D21g with transients in water content. Water Resour. Res. 51, 3303–3316.
Brusseau, M., Van Glubt, S., 2019. The influence of surfactant and solution composition
on PFAS adsorption at fluid-fluid interfaces. Water Res. 161 https://doi.org/
10.1016/j.watres.2019.05.095.
Brusseau, M.L., 2023. Determining air-water interfacial areas for the retention andtransport of PFAS and other interfacially active solutes in unsaturated porous media.
Sci. Total Environ. 884, 163730.
Brusseau, M.L., Guo, B., 2021. Air-water interfacial areas relevant for transport of per
and poly-fluoroalkyl substances. Water Res. 207, 117785.
Brusseau, M.L., Guo, B., Huang, D., Yan, N., Lyu, Y., 2021. Ideal versus nonidealtransport of PFAS in unsaturated porous media. Water Res. 202, 117405.
Carsel, R.F., Parrish, R.S., 1988. Developing joint probability distributions of soil water
retention characteristics. Water Resour. Res. 24, 755–769.
Feddes, R.A., Kowalik, P.J., Zaradny, H., 1978. Simulation of field water use and crop
yield. Simul. Monogr. PUDOC.
Felizeter, S., Juerling, H., Kotthoff, M., De Voogt, P., McLachlan, M., 2021. Uptake ofperfluorinated alkyl acids by crops: results from a field study. Environ Sci Process
Impacts 23. https://doi.org/10.1039/D1EM00166C.
van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulicconductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44, 892–898.
Gumuła-Kawęcka, A., Jaworska-Szulc, B., Szymkiewicz, A., Gorczewska-Langner, W.,
Pruszkowska-Caceres, M., Angulo-Jaramillo, R., Simunek, J., 2022. Estimation ofgroundwater recharge in a shallow sandy aquifer using unsaturated zone modelingand water table fluctuation method. J. Hydrol. 605, 127283 https://doi.org/
10.1016/j.jhydrol.2021.127283.
Guo, B., Zeng, J., Brusseau, M.L., 2020. A mathematical model for the release, transport,
and retention of per- and polyfluoroalkyl substances (PFAS) in the vadose zone.
Water Resour. Res. 56 (2), e2019WR026667.
Guo, B., Zeng, J., Brusseau, M.L., Zhang, Y., 2022. A screening model for quantifying
PFAS leaching in the vadose zone and mass discharge to groundwater. Adv. Water
Resour. 160, 104102.
Held, T., Reinhard, M., 2020. Project No. (FKZ) 3717 76 231 0. Report No. FB000332/
ENG. Remediation management for local and widespread PFAS contaminations.
Hutson, J.L., 2003. LEACHM: Leaching Estimation and Chemistry Model; a Process-Based
Model of Water and Solute Movement, Transformations, Plant Uptake and Chemical
Reactions in the Unsaturated Zone; Version 4. Cornell Univ, Ithaca.
Jackson, R.B., Canadell, J., Ehleringer, J.R., Mooney, H.A., Sala, O.E., Schulze, E.D.,
1996. A global analysis of root distributions for terrestrial biomes. Oecologia 108(3), 389–411.
Jobbágy, E.G., Jackson, R.B., 2000. The vertical distribution of soil organic carbon andits relation to climate and vegetation. Ecol. Appl. 10 (2), 423–436.
Johnson, G., 2022. PFAS in soil and groundwater following historical land application of
biosolids. Water Res. 211, 118035 https://doi.org/10.1016/j.watres.2021.118035.
Li, Y., Oliver, D.P., Kookana, R.S., 2018. A critical analysis of published data to discern
the role of soil and sediment properties in determining sorption of per andpolyfluoroalkyl substances (PFASs). Sci. Total Environ. 628, 110–120.
Minnesota Pollution Control Agency, 2024. Minnesota Department of Health, DNR &MDA. Authors Greene S. & Neuschler C. (Published February 2021). Minnesota's
PFAS Blueprint. Document number: p-gen1–22 Accessed on Jan 19th. https://www.
pca.state.mn.us/sites/default/files/p-gen1-22.pdf.
Mualem, Y., 1976. A new model for predicting the hydraulic conductivity of unsaturated
porous media. Water Resour. Res. 12 (3), 513–522.

11