Multidimensional simulation of PFAS transport and leaching in the vadose zone: impact of surfactant-induced flow and soil heterogeneities

PFAS are emergent contaminants of which fate and transport in the environment remain poorly understood. As surfactants, adsorption at air-water interfaces and solid surfaces in soils complicates the retention and leaching of PFAS in the vadose zone. Recent modeling studies accounting for the PFAS-specific nonlinear adsorption processes predicted that the majority of long-chain PFAS remain in the shallow vadose zone decades after contamination ceases—in agreement with many field measurements. However, some field investigations show that longchain PFAS have migrated to tens to a hundred meters below ground surface. These discrepancies may be attributed to model simplifications such as a one-dimensional (1D) representation of the homogeneous vadose zone. Another potentially critical process that has not been fully examined by the 1D models is how surfactant-induced flow (SIF) influences PFAS leaching in multidimensions. We develop a new three-dimensional model for PFAS transport in the subsurface to investigate the multidimensional effects of SIF and soil heterogeneities. Our simulations and analyses conclude that 1) SIF has a minimal impact on the long-term leaching of PFAS in the vadose zone, 2) preferential flow pathways generated by soil heterogeneities lead to early arrival and accelerated leaching of (especially long-chain) PFAS, 3) the acceleration of PFAS leaching in high water-content preferential pathways or perched water above capillary barriers is more prominent than conventional contaminants due to the destruction of air-water interfaces, and 4) soil heterogeneities are among primary sources of uncertainty for predicting PFAS leaching and retention in the vadose zone.


Introduction
Per-and polyfluoroalkyl substances (PFAS) are a group of synthetic chemicals that have been widely used since the 1940s in many industrial and consumer products, including nonstick and stain-resistant coatings, food containers, paper products, and fire fighting foams (Buck et al., 2011;ITRC, 2018;Wang et al., 2017). Resulting from large-scale manufacturing and wide use of PFAS in the past 80 years, they are now ubiquitous in the environment (e.g., soils, groundwater, sediments, and surface water) and have become emerging contaminants of international concern (e.g., Brusseau et al., 2020b;Hatton et al., 2018;Hu et al., 2016;Sharifan et al., 2021). The scientific understanding of the fate and transport of PFAS in the environment, especially in source zones under water unsaturated conditions, is still in its infancy. A growing body of field data demonstrates that PFAS concentrations in soils are orders of magnitude higher than those in groundwater beneath the source zones (up to seven orders of magnitude) and that the majority of PFAS remains near the top of the vadose zone (Anderson et al., 2016;Brusseau, et al., 2020b;Dauchy et al., 2019;McGuire et al., 2014;Xiao et al., 2015;Adamson et al., 2020). However, it remains poorly understood as to what primary mechanisms are controlling the long-term retention of PFAS in the vadose zone and the subsequent mass discharge to groundwater over time.
Many PFAS are surfactants (Kissa, 2001), thus they tend to adsorb at air-water interfaces and solid surfaces in soils and may stay in the vadose zone for long times before contaminating groundwater (Brusseau, 2018;Guo et al., 2020). Using surface tension data of two primary PFAS and measured air-water interfacial (AWI) areas in soil materials, Brusseau (2018) examined the impact of PFAS adsorption at air-water interfaces as a potential retention process in soils. More surface-tension-based theoretical analyses covering a wider range of PFAS and solution chemistry conditions have since been reported (Brusseau et al., 2019a(Brusseau et al., , 2019bBrusseau & Van Glubt, 2019;Costanza et al., 2019;Silva et al., 2019;Schaefer et al., 2019;Costanza et al., 2020). The impact of AWI adsorption has also been investigated by miscible-displacement experiments using packed columns of different soil media and water saturation conditions (Brusseau, et al., 2019a;Lyu & Brusseau, 2020;Lyu et al., 2018). The measured breakthrough curves were simulated by a onedimensional (1D) solute transport model, assuming steady-state flow under a uniform water saturation, that accounts for rate-limited and nonlinear adsorption at solid surfaces and air-water interfaces (Brusseau, 2020;Brusseau, et al., 2019a). The above surface-tension-based analyses, miscible-displacement experiments, and the 1D solute transport simulations all focused on laboratory conditions. Mathematical modeling of PFAS transport and retention under field conditions in the vadose zone has only begun recently. Guo et al., (2020) reported a new mathematical model that represents a comprehensive set of PFAS-specific transport processes including surfactant-induced flow (SIF), nonlinear and rate-limited solid-phase and AWI adsorption under transient variably saturated water flow. Numerical simulations of PFOS (Perfluorooctanesulfonic acid, C8HF17O3S) released at a model aqueous film-forming foam (AFFF)-impacted fire training area (FTA) site and the subsequent transport and retention in the vadose zone demonstrated that PFOS can be strongly retained in the vadose zone even under highly transient variably saturated flow driven by timedependent rainfall infiltration. Shortly after, a model similar to that of Guo et al. (2020) was also independently reported by Silva et al. (2020aSilva et al. ( & 2020b, where example simulations were presented in both 1D and 2D domains. The numerical simulations of PFAS leaching in the vadose zone reported by Guo et al. (2020) are generally consistent with soil PFAS concentration profiles measured at contaminated sites.
Namely, concentrations of long-chain PFAS in the soil profiles at many source zones show that the majority of the total mass stays in the shallow vadose zone (i.e., the first meter or two) Dauchy et al., 2019;Davis et al., 2007;Filipovic et al., 2015;Hale et al., 2017;Shin et al., 2012), and the concentration decreases exponentially along with the depth . However, despite that most of the long-chain PFAS were reported to stay in the shallow vadose zone, they are observed in deep vadose zones and groundwater at many sites even tens or a hundred meters below the land surface Dauchy et al., 2019;AFW, 2019). This discrepancy between the model simulations and the field data implies that additional processes and factors not incorporated in the prior modeling studies are likely influencing the leaching of PFAS in the vadose zone. A potentially important factor is the subsurface heterogeneities such as macropores, high-conductivity channels or fractures, soil lenses or layers, which were widely reported to provide preferential pathways for the transport and leaching of other solutes and contaminants in the vadose zone and groundwater (e.g., Brussea & Rao, 1990;Bekhit & Hassan, 2005;Köhne et al., 2006Köhne et al., , 2009Kung, 1990aKung, , 1990bŠimůnek et al., 2003;Zhou & Zhan, 2018).
Another important factor that has not been fully examined is the impact of SIF on PFAS leaching. Prior experimental and modeling studies on non-PFAS surfactants showed that surfactant-induced capillary pressure gradient can strongly influence variably saturated flow and solute transport in the vadose zone Gillham, 1994, 1999;Henry et al., 2002;Karagunduz et al., 2015). It is anticipated that similar effects of SIF may be observed for PFAS in source zones. Specifically, decreases in surface tension resulting from PFAS being present in porewater can cause induced drainage that on the one hand accelerates the transport of PFAS in the vertical direction due to increased porewater velocities, while on the other hand also leads to lateral spreading of PFAS that can potentially reduce vertical leaching. Concomitantly, surfactantinduced drainage also reduces water saturation resulting in an increase in the AWI area, which would then increase the retention of PFAS due to stronger AWI adsorption. It remains unclear as to how the above mechanisms will manifest to influence PFAS leaching in the vadose zone. Systematic investigations of the impact of both soil heterogeneity and SIF on PFAS leaching require more advanced models that represent flow and transport processes in multiple dimensions.
To address the critical questions regarding SIF and soil heterogeneity, we extend the model formulation of Guo et al. (2020) to develop a three-dimensional (3D) model that represents the transient variably saturated flow, SIF, and the comprehensive set of PFAS-specific transport processes in the vadose zone and groundwater. We apply the new model to simulate contamination scenarios at a model AFFF-impacted FTA site. A set of comprehensive numerical experiments are designed to examine the impact of SIF and soil heterogeneities on PFAS leaching and mass discharge to deep vadose zones or groundwater. Several types of typical soil heterogeneities are considered including soil lenses, high-conductivity pathways, and randomly distributed soil media in space. To the best of our knowledge, this is one of the first 3D models for PFAS transport in the vadose zone and multidimensional numerical simulation study of PFAS retention and leaching in the vadose zone accounting for a wide range of soil heterogeneities.
The multidimensional mathematical formulations and the numerical framework are presented in section 2. In section 3, data and parameters used for modeling PFAS transport are introduced. This is followed by section 4 that presents detailed descriptions of numerical experiments and analyses of the simulation results. In section 5, we discuss the implications of our findings. We close with concluding remarks in section 6. The variably saturated water flow in the subsurface may be described by the Richards' equation (Richardson, 1921;Richards, 1931), which incorporates the extended Darcy's law for variably saturated flow (Buckingham, 1907) in the mass balance equation of the water phase. Here, we present the Richards' where is the volumetric water content (cm 3 /cm 3 ). h is the water pressure head (cm). t is time (s). The vertical axis z is positive downward. is the tensor for saturated hydraulic conductivity and (-) is the relative permeability parameterized as a function of . The relative permeability function of Mualem (1976) and van Genuchten (1980) is used. We assume that the off-diagonal entries of are zero (Kij = 0 for i ≠ j) and that the porous medium is isotropic (K11 = K22 = K33 = Ks).

PFAS transport
The transport and retention of a PFAS in the subsurface may be described by the 3D advection-dispersion equation coupled with adsorption at the solid surfaces and air-water interfaces (e.g., Guo et al., 2020) where C is the aqueous concentration (µmol/cm 3 ). is the bulk density of the porous medium (g/cm 3 ). = / is the vector for the interstitial porewater velocity (cm/s), where = − (ℎ − ) is the Darcy flux (cm/s). The hydrodynamic dispersion tensor is the sum of the mechanical dispersion and molecular diffusion coefficients (cm 2 /s), which for an isotropic porous medium can be given by = | | + ( − ) | | + 0 (Bear, 1988), where is the identity matrix, (cm) and (cm) are longitudinal and transverse dispersivities, respectively; 0 is the molecular diffusion coefficient in free water (cm 2 /s), = 7/3 / 2 is the tortuosity factor for the water phase (Millington & Quirk, 1961).

Coupling between the water flow and PFAS transport
When PFAS is present in porewater, the variably saturated water flow (Eq. (1)) and PFAS transport (Eq. (2)) are fully coupled. Transient variably saturated flow drives advection and dispersion of PFAS and also influences the AWI adsorption by changing the AWI area. Concomitantly, changes in PFAS concentration modify the surface tension and capillary forces, and in turn influence the variably saturated water flow.
We parameterize as a function of ℎ following the soil-water characteristic (SWC) function of van Genuchten (1980). With PFAS in the porewater, we use the Leverett J function to obtain a new pressure head ℎ ′ = 0 • cos 0 cos ℎ that accounts for the change of surface tension (Leverett, 1941;Smith and Gillham, 1994;Henry & Smith, 2003;Guo et al., 2020). is the water contact angle, is the surface tension in the presence of PFAS in porewater, and ℎ 0 , 0 , and 0 are, respectively, the water pressure head, water contact angle, and surface tension when the PFAS concentration is zero. Replacing ℎ in the original SWC function with the new water pressure head ℎ ′ leads to a modified SWC function in the presence of PFAS in soil porewater (Guo et al., 2020).

Numerical framework
We present a numerical framework to solve the mathematical formulations presented in section 2.1 (i.e., Eqs. (1) and (2)) subject to initial and boundary conditions. We employ the cellcentered finite difference (in space) and backward Euler (in time) approximations to discretize both Eqs. (1) and (2). The mass conservative numerical scheme of Celia et al. (1990) is used for Eq. (1). For Eq. (2), a first-order upwind scheme is used for the advection term, while a central difference scheme is used for the dispersion term. The top boundary of the domain allows for surface evaporation and ponding conditions using an implementation similar to that reported in Guo et al. (2020).
We implement two numerical schemes for the coupling between Eqs.
(1) and (2) within each numerical time step. One employs a fully implicit (FI) framework that solves the two equations simultaneously, while the other uses a sequential implicit (SI) framework that iterates between the two equations until convergence. In either case, a Newton-Raphson solver is employed to solve the resulting nonlinear system of equations. The FI scheme assembles the discretized nonlinear formulations of Eqs. (1) and (2) and then linearizes by deriving a single Jacobian matrix to solve for ℎ and simultaneously at each iteration level. Conversely, the SI scheme linearizes the discretized nonlinear formulations of Eqs. (1) and (2) separately and generates two smaller Jacobian matrices-one for ℎ and the other for . An additional iteration loop is employed between the solvers for ℎ and at each numerical time step to ensure convergence. More details for the FI and SI frameworks are presented in Appendix A. We comment on the convergence properties of the two coupling frameworks in section 5 for the variably saturated flow and PFAS transport simulations examined in the present work. Our overall numerical implementation utilizes the gridding and automatic differentiation utilities of the Matlab Reservoir Simulation Toolbox (Lie, 2019). To validate the 3D implementation, we have compared our model to a widely-used open source simulator PFLOTRAN (Hammond et al., 2014). Excellent agreement between the two models is observed for both transient variably saturated flow and solute transport in 3D (see Figure  B.1 in Appendix B).

Data and parameters
We employ the newly developed 3D model to examine the impact of SIF and soil heterogeneity on PFAS leaching in the vadose zone. Three types of heterogeneities including the presence of soil lenses or layers, high-conductivity pathways, and randomly distributed soil media are examined. Four different porous media (covering a wide range of material properties) and two representative PFAS (i.e., one is short-chain and the other is long-chain) are used to construct the numerical experiments. In the following subsections, we introduce the detailed data and parameters for these porous media and PFAS.

Porous media
The four porous media that will be used for the numerical simulations are referred to as Accusand, Vinton, Hayhook, and high-K sand (i.e., a sand that has a greater hydraulic conductivity than the other three porous media). Detailed measured data and parameters are available for the first three from the literature (Guo et al., 2020;Jiang et al., 2020bJiang et al., , 2020aPeng & Brusseau, 2005). The last high-K sand is a porous medium that is hypothetically constructed to represent the preferential flow pathways. Accusand is a commercially available natural quartz sand (UNIMIN Corp.) with a median grain diameter of 0.35 mm and a uniformity coefficient of 1.1. The uniformity coefficient is defined as d60/d10 where d60 and d10 are the grain diameters that are greater than 60% and 10% of the sample by weight, respectively. Vinton and Hayhook were both collected locally in Tucson, Arizona. Vinton is a loamy sand. It has a median grain diameter of 0.23 mm and a uniformity coefficient of 2.4. Hayhook is a loam. Its mean grain diameter is 0.26 mm and the uniformity coefficient is 16. The difference between the uniformity coefficients indicates that Hayhook has a much wider range of particle size than that of Vinton. Even though Vinton and Hayhook have similar median grain sizes, Hayhook has many more smaller grains than Vinton. For example, the clay content of Hayhook (10%) is much higher than that of Vinton (3%).
The detailed parameters for the four porous media are summarized in Table 1. The capillary pressure as a function of water content, as well as the unsaturated hydraulic conductivity as a function of capillary pressure, are given in Figure C.1 in Appendix C. Most of the parameters for Accusand, Vinton, and Hayhook were determined in the laboratory. The remaining parameters whose measured values are not available are estimated as follows. Experimental observations and pore-scale models have demonstrated that the longitudinal and transverse dispersivities ( and ) are functions of the water content under unsaturated conditions (e.g., Bear & Cheng, 2010;Kirda et al., 1973;Maraqa et al., 1997;Padilla et al., 1999;Raoof & Hassanizadeh, 2013;Sato et al., 2003;De Smedt & Wierenga, 1984), but general functional forms have not become available. Here, we do not consider the dependence on water saturation and simply use the approximation for the longitudinal dispersivity under saturated conditions as = 83(log10(L/100)) 2.414 (Xu & Eckstein, 1995), where L is the apparent length scale. In the present study, we take L as the vertical depth of the vadose zone. The transverse dispersivity was reported to be generally smaller than (e.g., Bear & Cheng, 2010;Yule & Gardner, 1978).
≈ /10 is assumed in our simulations (Rockhold et al., 2016). Note that because the estimated and do not consider the impact of water saturation, they may lead to underestimation of dispersion in the vadose zone. The AWI area (Aaw) as a function of water saturation (Sw) for both Accusand and Vinton was fitted to measured values obtained by the aqueous-based interfacial partitioning tracer test (IPTT) and mass balance surfactant-tracer methods Guo et al., 2020), both of which were reported to be more representative than the AWI area measured by gas-based tracer methods for transport processes in the aqueous phase (Lyu et al., 2018;Brusseau et al., 2019;Guo et al., 2020). For Hayhook, only Aaw data obtained by the gas-based IPTT was available (Peng & Brusseau, 2005). The gas-IPTT-measured AWI areas are generally greater than those measured by the aqueous-IPTT method because a gas-phase surfactant tracer can access more AWI areas (likely due to its greater mobility) Jiang et al., 2020a). The measured aqueous-IPTT Aaw data for Accusand and Vinton are presented in Figure 1. The measured gas-IPTT Aaw data for Accusand, Vinton, and Hayhook are presented in Figure C.2 in Appendix C. Assuming that the Aaw measured by the gas-IPTT and the aqueous-IPTT differ by constant ratios for a porous medium, we estimate the aqueous-IPTT Aaw for Hayhook as Aaw,Hayhook = 2.3Aaw,Vinton. To test if the scaling estimate is reasonable, we then independently estimated the aqueous-IPTT Aaw for Accusand as Aaw,Accusand = 0.46Aaw,Vinton. The independently estimated aqueous-IPTT Aaw for Accusand was shown to match with the measured values well, which demonstrates that the scaling estimate appears to be a reasonable approach. When no directly measured data are available, other approaches may be used to determine AWI area, such as the thermodynamic approach that estimates the AWI area using the SWC (Morrow, 1970;Bradford & Leij, 1997;Bradford et al., 2014). However, we note that the AWI area estimated using the thermodynamic approach is generally lower than that measured by the interfacial tracer methods (Jiang et al. 2020), which is also the case for the soil media and the interfacial tracer measured AWI areas used in the present study.
The high-K sand is assigned a saturated hydraulic conductivity Ks that is significantly higher than that of the other three porous media, e.g., it is approximately 78.5 times the Ks for Accusand. The saturated water content and the residual water content are set to = 0.4 and = 0.02, respectively. The SWC curve and relative permeability are assumed to also follow the van Genutchen-Mualem models with = 0.1 (cm -1 ) and n = 5. The organic carbon content foc is set to one-tenth of that for Accusand. The Aaw for the high-K sand is supposed to be significantly lower than those of the other three porous media. We assume that its Aaw is one-tenth of that for Accusand at the same water saturation, leading to a maximum Aaw of 63.16 cm 2 /cm 3 at Sw = 0.002. The above parameters for the high-K sand are considered as the base case. To examine the sensitivity of PFAS leaching to the two parameters Ks and Aaw of the high-K sand, we consider a wider range of values for Ks (1/3Ks to 3Ks) and Aaw (3Aaw to 1/3Aaw) and compare them with the base case. Figure 1 The AWI area (Aaw) as a function of water saturation (Sw) for the four porous media (Accusand, Vinton, Hayhook, and high-K sand) used in the present study. The measured data points were obtained by aqueous-based tracer methods Guo et al., 2020). The coefficients for the fitted and estimated Aaw-Sw curves are presented in Table 1. 0.04 0.1 0.11 0.004 Note: All of the parameters were determined in the laboratory except for the ones denoted with superscript "*", which were estimated because measured values are not available.

Representative PFAS
We consider two representative PFAS-PFOS (long-chain) and PFPeA (Perfluoropentanoic acid, C5HF9O2, short-chain)-that have been widely found at contamination sites (Baduel et al., 2017;Brusseau, Anderson, et al., 2020;Dauchy et al., 2017;Høisaeter et al., 2019;Houtz et al., 2013;Sepulvado et al., 2011). As a long-chain PFAS, PFOS has a stronger adsorption capacity at both solid surfaces and AWI interfaces than that of PFPeA. The physicochemical and transport properties of the two PFAS including the Szyszkowski parameters for surface tension, Freundlich isotherm parameters for solid-phase adsorption, and the molecular diffusion coefficients are presented in Table 2.
The Szyszkowski parameters a (µmol/cm 3 ) and b (-) were obtained by fitting to surface tension data measured in synthetic groundwater (SGW) (Brusseau & Van Glubt, 2019). The surface tension of SGW with no dissolved PFAS is 0 = 71 dyn/cm. The Freundlich isotherms for PFOS in both Accusand and Vinton were obtained from prior sorption experiments reported in Brusseau (2020) and Brusseau et al. (2019a). Because no solid-phase adsorption data are available for PFPeA, we estimated the Kf and N parameters from the measured Freundlich isotherms of PFOA using the QSPR method developed by Brusseau (2019). Similarly, we used the same QSPR method to estimate the Freundlich isotherms for PFOS and PFPeA in Hayhook and high-K sand. Because PFOS and PFPeA are both anionic PFAS for which the solid-phase adsorption is dominated by hydrophobic interactions with the organic carbon content (Brusseau, 2018;Guelfo & Higgins, 2013), we assume that Kf scales with the organic carbon content among the different porous media. For example, for PFOS in Hayhook, ,Hayhook = ,Vinton ,Hayhook ,Vinton . The N values for Hayhook and High-K are assumed to be the same as those for Vinton and Accusand, respectively. The molecular diffusion coefficient D0 (cm 2 /s) for PFOS and PFPeA were obtained from laboratory measurements reported in Schaefer et al. (2019). Finally, because the change of contact angle is expected to be relatively small for anionic surfactants especially at relatively low aqueous concentrations relevant to PFAS-contaminated sites (Desai et al., 1992;Kissa, 2001;Rosen & Kunjappu, 2012), we assume that the contact angle does not change, i.e., ≈ 0 and cos / cos 0 ≈ 1. Note that for cationic, zwitterionic, or nonionic PFAS that interact more strongly with soil grain surfaces, the change of contact angle may need to be accounted for especially at higher aqueous concentrations (Desai et al., 1992;Henry & Smith, 2003). Mw (g/mol) is the molecular weight of the PFAS molecules. 1.2×10 -5 * has a unit of (μmol/g)/(μmol/cm 3 ) N , and N is dimensionless.

Climatic boundary conditions
We employ the precipitation and evapotranspiration data from a humid region represented by the site of Silas Little (Clark, 2004) in New Jersey, USA. The data are from 01/01/2005 to 12/31/2014 for 10 years and are available from the AmeriFlux database (https://ameriflux.lbl. gov). The 30-min resolution of the original data is smoothened to an 8-hour resolution. The 10-year data were repeated 7 times to generate a data set covering 80 years. The mean annual precipitation and potential evapotranspiration are 1,066 mm and 674 mm, respectively. All of the precipitation is assumed as rainfall for simplicity, i.e., we do not distinguish snow from rainfall. The evapotranspiration data is treated as the potential evapotranspiration to determine the surface evaporation fluxes. Additional details for the determination of the evaporation fluxes were reported in Guo et al. (2020).

Numerical simulation and results
We apply the 3D model presented in section 2 to investigate the lateral spreading and vertical leaching of PFAS in the vadose zone as impacted by SIF and soil heterogeneities. We consider the release of PFAS to the vadose zone at a model FTA site and construct two sets of scenarios--one involves a homogeneous vadose zone and the other considers several typical types of heterogeneities in a heterogeneous vadose zone.
The following is applied to all of the scenarios: (1) The FTA to which AFFF solutions are directly released has dimensions of 30 m × 30 m. (2) The vadose zone is assumed to have a depth of 4 m. The depth to groundwater table is assumed to be deeper than 4 m, so a free-drainage condition for water flow and a zero-gradient concentration condition for PFAS transport are employed at the bottom boundary. (3) Two periods in time are simulated including an activecontamination period of 30 years when the release of AFFF solutions occurs regularly due to fire training and a subsequent post-contamination period of 50 years during which the fire training stops and no PFAS are released to the vadose zone. (4) Real-world high-resolution rainfall and evapotranspiration data are applied at the top boundary (see more details of the data in section 3.3). During the active-contamination period, the AFFF solution is applied to the FTA once per 10 days representing a fire training event. The training frequency is within the range of approximately 12 to 44 times per year reported in the literature (e.g., Rotander et al., 2015;APEX, 2017). Each application is assumed to last for 30 minutes (at a rate of 22 mm/d). The AFFF solution released from fire training makes up 17 mm of water per unit area added to the FTA each year. The diluted AFFF solution is assumed to contain PFOS and PFPeA at concentrations of 100 mg/L and 0.528 mg/L (Høisaeter et al., 2019), respectively, if not otherwise specified. (5) The initial water pressure head of the domain is h0 = -100 cm and the domain is assumed initially PFAS-free, i.e., C0 = 0. (6) Convergence tests have been conducted for the numerical discretizations and all of the numerical solutions have converged.
In the following subsections, we present more specifics for the scenarios, which will be followed by detailed simulation results and analysis.

Scenario 1: Impact of surfactant-induced flow on PFAS leaching
In this scenario, we focus on examining the impact of SIF on PFAS leaching in a homogenous vadose zone. In particular, we quantify 1) the potential acceleration of vertical leaching and 2) the strength of lateral spreading resulting from the SIF. As shown in Figure 2, we consider a crosssectional domain (30 m × 4 m) of a homogeneous vadose zone represented by one of the three porous media (i.e., Accusand, Vinton, and Hayhook). We assume that the FTA is symmetric and thus only simulate half of the FTA by imposing a non-flux condition on the left boundary. Firetraining occurs at the upper boundary covering 0 ≤ x ≤ 15 m. In addition to considering PFOS of 100 mg/L and PFPeA of 0.528 mg/L in the AFFF solution, we also simulate another case with PFOS at a much higher concentration of 1,000 mg/L, which can be considered as an upper limit for the PFOS concentration in the diluted AFFF solutions applied to FTAs (Houtz et al., 2013;Guo et al., 2020). Finally, a non-reactive tracer at a concentration of 0.528 mg/L is also simulated for comparison. The 2D domain is discretized using uniform grid cells ∆x = ∆z = 10 cm. We present the 2D PFAS plumes at t = 80 years (Figure 3). Because the presence of (longchain) PFOS in porewater leads to much stronger SIF than that of (short-chain) PFPeA, we only present the plumes for PFOS. The simulation results from both the release concentrations of 100 mg/L and 1,000 mg/L in all three porous media are presented. The lateral spreading of PFOS is shown to be generally small. The extent of spreading in the three porous media has the following order of Accusand > Vinton > Hayhook and it is greater under the higher release concentration of 1,000 mg/L. The largest extent of spreading occurs when 1,000 mg/L of PFOS is applied to the Accusand--the plume extends 1.2 m beyond the vadose zone directly underneath the FTA (x > 15 m), which is still insignificant compared to the overall size of the plume in the lateral direction. To quantify the impact of SIF, we turn off the change in surface tension (i.e., set = 0 in Equation (3)) and perform another set of simulations for comparison. The results show that SIF has a relatively minor impact, i.e., the PFOS plumes with and without accounting for SIF are similar ( Figure 3). To further quantify the impact of lateral spreading on PFAS leaching, we also conduct 1D simulations with and without accounting for SIF, similar to those reported in Guo et al., 2020, and compare the vertical concentration profiles with our 2D simulations ( Figure 4; the concentrations from the 2D simulations are averaged along the horizontal dimension). The difference between the 2D and 1D simulations is almost negligible, confirming that the lateral spreading due to SIF has a minimal impact on PFAS leaching.
In addition to SIF, molecular diffusion and transverse mechanical dispersion can also cause lateral spreading. To evaluate the extent of lateral spreading contributed solely by SIF, we compute and compare the mass fraction of PFAS that migrates outside of the vadose zone directly underneath the FTA with and without accounting for SIF for both PFOS and PFPeA ( Figure 5). The non-reactive tracer is also presented as a reference. We make the following observations: 1) The lateral spreading is mainly driven by transverse mechanical dispersion (due to strong porewater velocities in the vertical direction), rather than SIF. Note that molecular diffusion is negligible compared to the transverse mechanical dispersion. 2) Consistent with Figure 3, the extent of lateral spreading follows the order of Accusand > Vinton > Hayhook. This is because the lateral spreading is mainly caused by transverse mechanical dispersion and that the vertical porewater velocities for Accusand > Vinton > Hayhook (see their Ks in Table 1). 3) Among the three porous media, the lateral spreading due to SIF is the weakest in Hayhook. A closer inspection reveals that the aqueous concentration of PFOS is the smallest in Hayhook-the majority of PFOS is adsorbed at solid surfaces and air-water interfaces due to strong solid-phase and AWI adsorption. A lower aqueous concentration leads to weaker SIF. For example, the maximum aqueous concentration is 11 mg/L for Accusand, 10 mg/L for Vinton, and 4.5 mg/L for Hayhook at t = 10 years when PFOS is applied at 1,000 mg/L. 4) Though lateral spreading due to SIF is measurable for the release concentration 1,000 mg/L for PFOS (approximately 5.5%), it is much smaller for PFOS at the release concentration of 100 mg/L (< 1.2% for all three porous media) and also for PFPeA (< 0.1% for all three porous media). In particular, the spreading of PFPeA is almost identical to that of the non-reactive tracer, indicating that SIF has no impact on the lateral spreading of PFPeA.   In this scenario, we simulate PFAS leaching in the vadose zone in the presence of three typical types of soil heterogeneities. To accurately represent the heterogeneities in our numerical discretization, we focus on subdomains of a vadose zone directly underneath the FTA to allow for a higher spatial numerical resolution. The schematics and the spatial arrangement of the heterogeneity inclusions for the three test cases are presented in Figure 6. More specific information about the heterogeneities for each set of the test cases is provided below.

Soil layering and lens
A soil lens or a fully extended layer is embedded in an otherwise homogeneous vadose zone composed of the Vinton soil (Figure 6a). An x-z cross-sectional domain of 5 m × 4 m is considered. The soil lens is elliptic (center: (x, z) = (2.5 m, 1.5 m), major axis: 4 m; minor axis: 1 m). Alternatively, we consider a fully extended horizontal layer at 1 ≤ z ≤ 2 m. The lens or layer is represented by one of the three porous media, i.e., Accusand, Vinton, and Hayhook, leading to a set of five test cases: homogeneous Vinton, Vinton with an embedded lens or layer of Accusand, and Vinton with an embedded lens or layer of Hayhook. The 2D domain is discretized using uniform grid cells ∆x = ∆z = 5 cm.
We present the time-dependent vertical total concentration profiles ( Figure 7) and the cumulative breakthrough curves (CBTCs) (Figure 8), i.e., the cumulative mass of PFOS and PFPeA leaching through the bottom of the domain, for all five test cases. The 1D vertical concentration profiles were obtained by averaging the 2D plumes along the x-axis. In addition, we also present the spatial distributions of water pressure head (h), water saturation (Sw), AWI area (Aaw), and total concentration (Ctot) (Figure 9).
Comparisons of the vertical concentration profiles in Figure 7 show that the inclusion of a heterogeneous soil lens or a layer in an otherwise homogeneous vadose zone significantly influences the leaching behavior of PFAS. This is partly due to the fact that the strengths of both solid-phase and AWI adsorption for PFAS in Vinton are stronger than that in Accusand, while weaker than that in Hayhook. Notably, the leaching of both PFOS and PFPeA is shown to be faster in the presence of a lens than that in the presence of a fully extended layer despite that they are composed of the same soil medium. These observations are confirmed by the arrival times (ta, defined as the time at which PFAS reaches the bottom of the domain) shown in the CBTCs ( Figure  8). For PFOS in the homogeneous Vinton, ta = 59 years. The presence of an Accusand lens leads to a significantly earlier arrival (ta = 34 years) that is even 19 years earlier than that for a fully extended Accusand layer (ta = 53 years). The presence of a lens and a fully extended layer of Hayhook leads to an arrival time of ta = 54 years and ta = 62 years, respectively. Even though the presence of either a Hayhook lens or a fully-extended layer enhances the overall retention ( Figure  8), it is rather interesting that a Hayhook lens leads to an arrival time of 5 years earlier (ta = 54 years) than that for the homogeneous case, which indicates the presence of some type of preferential transport. Similar trends discussed above can be observed in the arrival times of PFPeA. A closer inspection reveals that the above-mentioned different behaviors for the transport of PFAS between the inclusion of a lens and a fully-extended layer are caused by the so-called "funneled" flow (i.e.., one type of preferential flow) generated by the presence of the lens (Kung 1990a,b). For PFAS in particular, it is quite interesting that the funneled flow not only introduces preferential transport of PFAS by advection and dispersion but also amplifies leaching by reducing the strength of AWI adsorption. A more detailed analysis is provided below. Figure 7 Time-dependent profiles of the total concentrations (Ctot) for PFOS and PFPeA as influenced by the presence of a soil lens or a fully extended layer. The Ctot of PFOS or PFPeA is defined as its total mass (including the mass in the aqueous phase and that adsorbed at air-water interfaces and solid surfaces) per unit volume of the porous medium. These 1D profiles are obtained by averaging the 2D plumes along the horizontal dimension. Note that the results for PFPeA are only presented up to t = 33 years when the majority of the mass had left the domain, and this applies to other figures for PFPeA in the rest of the paper. When a lens of Accusand is embedded in an otherwise homogeneous vadose zone comprised of Vinton, water accumulates above the Vinton-Accusand interface before it can enter the underlying Accusand (the matrix potential in the Accusand is greater than that in the Vinton; see the second plot in row 1 of Figure 9). The accumulated water then moves around the Accusand lens (partly facilitated by the slopes), which eventually leads to preferential flow pathways (see column 2 of Figure 9). A Hayhook lens also introduces funneled flow though the mechanism is slightly different-water moves around the Hayhook lens resulting from its low hydraulic conductivity that makes water difficult to infiltrate. Note that the preferential flow generated by the Hayhook lens is much less profound than that by the Accusand lens. Figure 9 Spatial distributions of water pressure head (h), water saturation (Sw), AWI area (Aaw), and total PFOS concentration (Ctot) at t = 78 years. Each column corresponds to a test case that involves a different soil heterogeneity.
Like other solutes, the funneled flow leads to preferential transport of PFAS by advection and dispersion. However, what is unique about PFAS is that the preferential flow leads to greater water saturation that subsequently destructs AWI areas along the pathway (see the second and fourth plots of row 3, Figure 9). Because the AWI area is a primary retention process for PFAS, especially the long-chain ones, the destruction of AWI areas on the preferential flow pathway further contributes to accelerating the leaching of PFAS. This mechanism can be illustrated by comparing the accelerated leaching of PFOS to that of PFPeA. Namely, AWI adsorption plays a much more important role in the retention of the long-chain PFOS, the funneled flow caused by the lens of Accusand and Hayhook leads to greater acceleration of leaching for PFOS than that for PFPeA as demonstrated by the arrival times presented in Figure 8.
The other interesting finding unique to PFAS is that the presence of Accusand and Hayhook (either as a lens or a fully extended layer) accelerates the leaching of PFOS in the region above them, but not for PFPeA. Figure 7 shows that the mass of PFOS in the first meter is smaller when a lens or a fully extended layer of Accusand or Hayhook is present. Conversely, the mass of PFPeA in the first meter increases in the presence of Accusand or Hayhook. This is because both Accusand and Hayhook act as barriers that lead to the accumulation of water above the Vinton-Accusand or Vinton-Hayhook interface (i.e., the water saturation is greater than that for the homogeneous case). The greater water saturation reduces the amount of Aaw, and as a result, decreases the retention of PFOS in the Vinton above. Conversely, because AWI adsorption does not play a major role in the retention of PFPeA, water accumulating at the Vinton-Accusand or Vinton-Hayhook interface leads to more PFPeA (in the aqueous phase) accumulating above Accusand or Hayhook.

High-conductivity pathways
We simulate PFAS leaching in the presence of well-connected high-conductivity pathways that are composed of the high-K sand (Figure 6b). All of these high-K sand channels have a width of 5 cm. In addition to the base simulations using the parameters presented in Table 1, we also perform a sensitivity analysis by conducting additional simulations using a wider range of values for Ks and Aaw of the high-K sand. 1/3 and 3 times of the base values of these two parameters are used in the sensitivity analysis. The domain is discretized using uniform grid cells ∆x = ∆z = 2.5 cm.
The time-dependent vertical total concentration profiles (averaged from the 2D plumes along the x-axis; Figure 10) and the CBTCs (Figure 11) demonstrate that the high-conductivity pathways lead to much earlier arrival of PFAS at the bottom of the domain, especially for the long-chain PFOS. PFOS arrives at the bottom of the domain at ta = 4 years in the presence of the highconductivity pathways, which is approximately 1/15 of that for the homogenous case (ta = 59 years). Due to minimal retention from solid-phase and AWI adsorption, the arrival time for the short-chain PFPeA is short for the homogeneous case (ta = 1.23 years). However, the presence of the high-conductivity pathways further reduces it to ta = 0.05 years.
Though the high-conductivity pathways lead to the early arrival of PFAS at the bottom of the domain, they also provide "shortcut circuiting" channels that can make the leaching of PFAS in the surrounding lower conductivity regions less effective, especially during the post-contamination period. This phenomenon can be observed in the vertical concentration profiles in Figure 10 as well as the 2D concentration distributions in Figure 12. The vertical concentration profiles show that more PFOS and PFPeA remain in the shallower vadose zone (0 ≤ z ≤ 1 m and 0 ≤ z ≤ 2.5 m, respectively) during the post contamination period (t > 30 years) in the presence of the highconductivity channels when compared to the homogeneous case. A closer inspection of the 2D concentration distributions (Figure 12) reveals that, during the active-contamination period, some PFAS have entered the lower conductivity regions driven by advection. Due to the presence of the preferential pathways provided by the high-K sand channels, leaching of PFAS in those regions is less effective compared to that of the homogenous case.
Finally, the sensitivity test cases using a wider range of values for the hydraulic conductivity and AWI area for the high-K sand show a minimal difference from the base simulations ( Figures  10 and 11). This demonstrates that the overall leaching behaviors for both PFOS and PFPeA observed above are robust for a relatively wide range of high-K sand provided that the high-K sand has a greater conductivity and weaker adsorption capacities than that of Vinton. Figure 10 Time-dependent total concentration profiles of PFOS and PFPeA with vs. without the presence of high-conductivity pathways. The Ctot of PFOS or PFPeA is defined as its total mass (including the mass in the aqueous phase and that adsorbed at air-water interfaces and solid surfaces) per unit volume of the porous medium. These 1D profiles are obtained by averaging the 2D plumes along the horizontal dimension. The gray lines indicate the concentration profiles corresponding to the sensitivity test cases. Figure 11 Cumulative breakthrough curves for PFOS and PFPeA over time with vs. without the presence of high-conductivity pathways. The gray lines indicate the concentration profiles corresponding to the sensitivity test cases. The cumulative mass is normalized by the total PFOS or PFPeA released to the vadose zone during the entire active-contamination period. Figure 12 2D plumes of the total concentrations for PFOS and PFPeA over time with vs. without the presence of high-conductivity pathways.

Random heterogeneity
In this final set of test cases, we simulate PFAS leaching in the presence of much more complex 3D heterogeneities that follow a random distribution (Figure 6c). The 3D domain has dimensions of 1 m × 1 m × 4 m and is discretized using uniform grid cells ∆x = ∆y = ∆z = 10 cm. We use the Gaussian Kriging variogram model in the gstat package (Pebesma, 2004) to generate 3D random fields with the following parameters: nugget = 0, sill = 10 cm, correlation lengths = 40 cm and 4 cm in the horizontal and vertical directions (the correlation lengths are within the range reported in the literature (e.g., Russo, 1991;Tseng and Jury, 1993)), the maximal number of nearest observations = 100, and maximal distance from the prediction location used for prediction = 50 cm. Note that the correlation lengths are used for illustrative purposes. Other correlation lengths will likely lead to different quantitative leaching behaviors, but as long as preferential flow pathways are present, accelerated PFAS leaching due to reduced air-water interfacial adsorption as demonstrated in sections 4.2.1 and 4.2.2 will occur. The random values are then normalized into the range from 0 to 1. We then uniformly divide the cumulative density function (CDF) into four intervals assign them to the four soils, i.e., High-K sand, Accusand, Vinton, and Hayhook, as demonstrated in Figure 6(c). Three different sets of fractions among the four porous media are considered-0%, 3%, and 10% are assigned to the high-K sand, respectively, and the remaining portions are equally assigned to the other three porous media. For each set of fractions, we generate 8 realizations to conduct sensitivity analysis. Finally, we note that because the spatial heterogeneities have been explicitly represented in the 3D domain, we employ a lower dispersivity coefficient ( = 2 cm) in each homogeneous soil medium in the simulations.
We present the time-dependent vertical total concentrations profiles (averaged from the 3D plumes along the x-y plane; Figure 13) and the CBTCs (Figure 14) for the 24 simulations conducted using the randomly generated heterogeneous domains that contain 0%, 3%, and 10% of high-K sand (each has 8 realizations). The simulation results conducted using the homogeneous porous media (i.e., Accusand, Vinton, and Hayhook) are also presented for comparison. Figure 13 Time-dependent total concentration profiles for PFOS and PFPeA in homogeneous vs. randomly heterogeneous vadose zones that contain different fractions of high-K sand (0%, 3%, and 10%). The total concentration of PFOS or PFPeA is defined as its total mass (including the mass in the aqueous phase and that adsorbed at air-water interfaces and solid surfaces) per unit volume of the porous medium. These profiles are obtained by averaging the 3D plumes along the x-y plane.
For PFOS, the presence of high-K sand has a significant impact on its long-term leaching behavior. With 0% of high-K sand, the differences among the concentration profiles from the 8 heterogeneous realizations are minor and all of them are bounded by the concentration profiles of the homogenous Hayhook and Accusand cases (Figure 13). The arrival times for the 8 realizations vary from 43 to 76 years. Only 0.1% to 1.9% of the total PFOS has been discharged from the bottom of the domain at t = 80 years ( Figure 14). However, the addition of 3% high-K sand greatly accelerates the leaching process and leads to much more significant variations among the 8 realizations. The arrival times are reduced to 10 to 66 years. The earliest arrival time is more than 24 years earlier than that for the homogenous Accusand case. For this realization, more than 93% of total PFOS has been discharged from the bottom of the domain at t = 80 years ( Figure 14). When the fraction of the high-K sand is further increased to 10%, the leaching becomes even faster, and much greater variations are observed among the 8 realizations. The arrival times are reduced to 3 to 27 years. In the realization where leaching occurs the fastest, 100% of the total PFOS has been discharged from the bottom of the domain at t = 40 years. These synthetic simulations illustrate that even the presence of a small fraction of high-conductivity media can greatly accelerate the leaching of PFOS (i.e., long-chain PFAS) in the vadose zone. Figure 14 Cumulative breakthrough curves for PFOS and PFPeA in homogeneous vs. randomly heterogeneous vadose zones that contain different fractions of high-K sand (0%, 3%, and 10%). The cumulative mass is normalized by the total PFOS or PFPeA released to the vadose zone during the entire active-contamination period.
Conversely, the presence of the high-K sand shows a minimal impact on the long-term leaching of the short-chain PFPeA, which is consistent with the findings reported in sections 4.2.1 and 4.2.2. This is because the short-chain PFPeA is minimally retained in any of the four porous media (Accusand, Vinton, Hayhook, and high-K sand), and as a result, the different heterogeneity realizations generated using these four porous media do not significantly change the overall retention capacity for PFPeA in the vadose zone.
Finally, we present the 3D plumes to further illustrate and visualize the impact of the high-K sand on the leaching of PFOS and PFPeA in the vadose zone. A realization with no inclusion of high-K sand and a realization with 10% of high-K sand (the one that leads to the fastest leaching) are shown as examples ( Figure 15). The comparisons between the 3D plumes confirm that the presence of high-K sand significantly accelerates the migration and leaching of the long-chain PFOS, while it has a minimal impact on that of the short-chain PFPeA. In the realization with 10% of high-K sand, the high-K sand forms well-connected channels that provide preferential pathways for the leaching of PFOS. However, these preferential pathways appear to be much less effective for PFPeA, primarily because the retention capacity of the other three porous media for PFPeA is already minimal as was discussed above. Figure 15 Comparison of the 3D plumes of the total concentrations for PFOS and PFPeA in two realizations of randomly heterogeneous vadose zones that contain 0% vs. 10% high-K sand.

Discussion
We have developed a new 3D model that represents the coupled variably saturated flow and transport processes relevant to PFAS retention and leaching in the vadose zone, including SIF and nonlinear adsorption of PFAS at the solid surfaces and air-water interfaces. The newly developed multidimensional modeling framework allowed us to conduct a series of numerical experiments to investigate the impact of two critical factors-SIF and soil heterogeneities-on the long-term leaching of PFAS in the vadose zone.
A prior study employing a 1D model (Guo et al., 2020) showed that the enhanced vertical drainage resulting from SIF has a negligible influence on the long-term leaching of PFAS. However, it remains unclear as to whether the lateral spreading of PFAS plumes caused by SIF will lead to potentially enhanced retention of PFAS in the vadose zone. In the present work, we demonstrate that SIF has a minimal impact, either via enhanced vertical drainage or lateral spreading, on the long-term leaching of PFAS. Our 2D numerical experiments (section 4.1) show that, even at the higher-end release concentration of 1,000 mg/L, only 5.5% of PFOS migrated outside the vadose zone directly underneath the FTA at t = 80 years. Notably, SIF only contributed to less than 1.3% of the 5.5%-the remaining 4.2% was caused by transverse mechanical dispersion in the lateral dimension. The reason for the minimal impact of SIF is the relatively low aqueous concentration of PFAS in porewater. Though PFAS can be released at relatively high concentrations especially at FTAs, after entering the vadose zone the majority of long-chain PFAS (the ones that are more surface active and can modify surface tension) are adsorbed at air-water interfaces and solid surfaces. This mass partitioning due to solid-phase and AWI adsorption significantly reduces the aqueous concentrations in porewater and as a result, reduces the change in surface tension. For example, the maximum aqueous concentrations of PFOS for the simulations conducted under all conditions in section 4.1 are less than 15 mg/L even though PFOS was released at a concentration of 1,000 mg/L. The minimal impact of the SIF has been further confirmed by excellent agreement among the vertical concentration profiles simulated by the 2D model and those by the 1D model with and without accounting for SIF.
The 2D and 3D numerical experiments for heterogeneous vadose zones (section 4.2) show that soil heterogeneity can be a primary factor that contributes to the early arrival of long-chain PFAS in groundwater. We illustrate that the presence of soil lenses and high-conductivity channels generates preferential pathways for flow and transport that can lead to the earlier arrival of PFAS compared to that of PFAS leaching in a homogeneous vadose zone. This suggests that the presence of soil heterogeneity is a plausible explanation for the early arrival of PFAS in the deep vadose zone or groundwater at tens to a hundred meters observed at some contaminated sites Dauchy et al., 2019;AFW, 2019). For example, Dauchy et al. (2019) reported deep seepage of both short-chain and long-chain PFAS to the groundwater below a depth of 15 m at an FTA site in France despite the presence of multiple clay layers. Similarly, at an Air Force Base in Tucson, AZ, long-chain PFAS reached a depth of over 70 m through a thick vadose zone (AFW, 2019). While the above field evidence shows a significant deviation of PFAS leaching rate from the simulation results based on homogeneous soils, the early arrivals appear to be consistent with our model simulations that have explicitly represented soil heterogeneities.
While preferential flow and transport resulting from soil lenses or high-conductivity pathways have been widely reported for other non-PFAS contaminants in the vadose zone and groundwater (e.g., Brusseau & Rao, 1994;Bekhit & Hassan, 2005;Köhne et al., 2006Köhne et al., , 2009Kung, 1990aKung, , 1990bŠimůnek et al., 2003;Zhou & Zhan, 2018), our simulations have revealed some interesting behaviors unique to PFAS. For example, the preferential flow pathways generated by the presence of the soil lens (i.e., "funneled" flow) not only leads to preferential transport of PFAS by advection and dispersion (like for nonreactive solutes) but also reduces the strength of retention for PFAS that further enhances preferential transport. The reduced retention, as analyzed in section 4.2.1, is caused by the increased water saturation along the preferential flow pathways that destructs AWI areas and hence decreases AWI adsorption. Another interesting behavior is that the spatial heterogeneity influences the long-term leaching of long-chain PFAS much more strongly than their short-chain counterparts, which is illustrated by the 2D experiments in sections 4.2.1-4.2.2 as well as the 3D numerical experiments conducted using randomly distributed heterogeneities in section 4.2.3. This difference between the long-chain and short-chain PFAS may be attributed to two factors, and both of which are a direct result of the long-chain PFAS being more strongly adsorbed at solid surfaces and air-water interfaces. Namely, for the long-chain PFAS, the preferential transport pathways are amplified (due to reductions in AWI area), and the differences in the adsorption capacities among the different soil media have a more significant impact on the overall retention. These PFAS-specific responses to soil heterogeneities imply that, compared to other more traditional contaminants (e.g., nonreactive solutes), more significant errors (especially for the long-chain compounds) may be introduced in the predictions of long-term leaching of PFAS if soil heterogeneities are not fully represented in the numerical model.
We comment on some assumptions employed in our simulations and how they may potentially influence the results. As presented in section 2.1.2, our model assumes equilibrium adsorption on solid grain surfaces and air-water interfaces. Prior miscible-displacement experiments indicate that rate-limited adsorption is likely insignificant under field-scale conditions where porewater velocities are much lower than the typical porewater velocities used in column experiments (Brusseau, 2020;Brusseau et al., 2019b). A recent study of batch experiments showed that kinetic solid-phase desorption may be more pronounced in field-aged soils, especially for shorter-chain PFAS (Schaefer et al., 2021). When rate-limited adsorption is present, PFAS leaching may be enhanced during the active contamination period (due to rate-limited adsorption) while reduced during the post-contamination period (due to rate-limited desorption).
The relative importance of solid-phase adsorption will depend on the specific soil media and PFAS under consideration. The foc values used in our study are relatively low compared to those for surface soils. However, since foc generally decreases rapidly along depth (Elzein & Balesdent, 1995;Jobbágy & Jackson, 2000), the foc especially those for Vinton soil and Hayhook soil should be representative for soils and sediments in the deeper vadose zone (approximately 1 m below the land surface). We recognize that foc in the first meter below the land surface may be much greater than the foc employed in the present study. Additionally, more complex sorption processes such as electrostatic interactions may lead to enhanced solid-phase adsorption in particular for cationic and zwitterionic PFAS (Barzen-Hanson et al., 2017;Mejia-Avendaño et al., 2020;Nickerson et al., 2020;Xiao et al., 2019). Under these conditions, our simulations would overestimate PFAS leaching.
The simulations conducted in the present study assume that each homogeneous soil medium is isotropic. Anisotropy with greater Kxx or Kyy than Kzz would in general lead to greater lateral spreading of the plume. In a homogeneous vadose zone, greater lateral spreading will reduce PFAS leaching in the vertical direction. While in a heterogeneous vadose zone, greater lateral spreading may help the PFAS find preferential flow pathways and hence enhance leaching. However, when heterogeneities are explicitly represented (like in the simulations presented in section 4.3), the domain becomes effectively anisotropic and hence would partially account for anisotropy. Another assumption employed in the simulations is the spatially uniform PFAS loading at the upper boundary. Spatial variation in PFAS loading rate may lead to spatially heterogeneous leaching rates. Locations with greater PFAS loading rates will likely have greater leaching rates resulting from reduced solid-phase and air-water interfacial adsorption.
Finally, we comment on the two numerical schemes (i.e., FI and SI) employed to couple the transient variably saturated flow and PFAS transport equations in our modeling framework. Both of these two coupling schemes are shown to provide convergent solutions. The nonlinear iteration loops for the FI method were shown to converge faster than that of the SI method for the simulations conducted in the present work, though the single linear Jacobian matrix in the FI scheme requires a preconditioner for improved convergence (if iterative linear solvers are used) and accuracy properties (Benzi, 2002). Our observations are consistent with a recent numerical analysis of the iterative schemes for the coupling between surfactant transport and variably saturated flow (Illiano et al. 2020). Note that AWI adsorption and solid-phase adsorption were not accounted for in the model formulation of Illiano et al. (2020).

Conclusions
In the present study, we have developed a multidimensional numerical modeling framework to simulate PFAS transport and leaching in the vadose zone. The model formulation represents transient variably saturated flow, PFAS advective and dispersive transport, nonlinear solid-phase and AWI adsorption, and SIF. The coupled flow and transport equations are solved using a fully implicit or a sequential implicit framework, either of which was shown to give convergent solutions. Employing the newly developed multidimensional modeling framework, we have conducted a series of numerical experiments at model FTA sites to understand and quantify the impact of SIF and soil heterogeneities on the long-term leaching of PFAS in the vadose zone. A long-chain (PFOS) and a short-chain (PFPeA) PFAS that are commonly found at contaminated sites are used as examples to examine the influence of chain length on the leaching behavior. All of the simulations last for 80 years including 30-year active contamination due to regular fire training activities and 50-year post contamination when fire training activities have been stopped. Our simulation results and analyses lead to the following conclusions: (1) SIF, resulting from the change of surface tension due to the presence of PFAS in porewater, has a minimal impact on the long-term retention and leaching of PFAS in the vadose zone. This has been confirmed by both 1D simulations that account for the induced drainage in the vertical dimension as well as 2D simulations that further represent the induced lateral spreading of the PFAS plume. The primary reason is that though PFAS may be released at high concentrations (e.g., at FTA sites), their aqueous concentration decreases significantly after entering the vadose zone due to strong adsorption at the solid surfaces and air-water interfaces. The relatively low aqueous concentrations have a relatively minimal impact on the porewater surface tension. Overall, our simulations imply that it may be acceptable to decouple the Richards' equation and the PFAS transport equations (i.e., assuming no SIF) when focusing on the long-term retention and leaching of PFAS in the vadose zone.
(2) Soil heterogeneity is a primary factor that controls the retention and leaching of PFAS in the vadose zone. The presence of heterogeneity generates preferential pathways for flow and transport that can lead to much earlier arrival of PFAS at a given depth compared to that in the homogenous vadose zones.
(3) A unique property differentiating PFAS from conventional non-PFAS contaminants is that preferential flow pathways can further amplify the preferential transport due to the destruction of air-water interfaces by the increased water saturation. We have also demonstrated that soil heterogeneities influence the long-term leaching of the long-chain PFAS more strongly than that of their short-chain counterparts.
(4) Given the significant impact of heterogeneity on PFAS leaching in the vadose zone, our study suggests that it is essential to employ multidimensional models for studying the PFAS retention and leaching at contaminated sites where complex heterogeneities are present.

A. The FI and SI frameworks
We introduce the numerical procedures for the FI and SI frameworks employed to solve the coupled equations in the present work. Let ( , ) = 0 and ( , ) = 0 represent the spatially and temporally discretized systems of Eqs. (1) and (2), respectively.

B. Numerical code validation
We validate our numerical implementation by comparing it to a widely-used open source simulator PFLOTRAN (Hammond et al., 2014). In the numerical experiment, at t = 0, a nonreactive tracer at a concentration of C0 = 0.5 mg/L is uniformly distributed in a 1 m × 1m × 1m 3D domain composed of Accusand. The domain is discretized using uniform cells ∆x = ∆y = ∆z = 1 cm. The water table depth is 90 cm and is kept constant during the simulation. The initial water pressure head is hydrostatic. The four side boundaries are all set as non-flux boundaries for water flow and solute transport. At the center of the top boundary, starting at t = 0, constant infiltration of clean water at 1 cm/d is applied to a 10 cm × 10 cm square area and lasts for 48 hours. The rest of the top boundary is set to zero flux. The 3D spatial distributions of water pressure head (h), water saturation (Sw), and aqueous concentration (C) are computed. The cross-sections of these distributions at y = 49.5 cm at the end of simulation using results computed from our code and PFLOTRAN are presented in Figure B.1. Excellent agreement between the results from the two codes is observed. Figure B.1 Comparisons of 2D cross-sections (at y = 49.5 cm) for water pressure head (h), water saturation (Sw), and aqueous concentration (C) between the results computed from our code and that from PFLOTRAN at the end of the simulation (t = 48 h).
C. Additional data and parameters for the four porous media Figure C.1 (a) Measured and fitted capillary pressure (-h) as a function of soil water content (θ) for the four porous media used in this study. Note that the data points near the oven-dry water content were removed to better fit the van Genuchten model for the range of water content that is more relevant to the field conditions. (b) Unsaturated hydraulic conductivity (K) as a function of capillary pressure head (-h).
Figure C.2 Measured and fitted AWI area (Aaw) as a function of water saturation (Sw = θ/θs). The measured data from aqueous-based tracer experiments are fitted to a second-degree polynomial function of water saturation (i.e., = 2 2 + 1 + 0 , see solid lines). The data from gas-IPTT experiments are fitted using the function and parameters provided by Peng & Brusseau (2012) (see dotted lines).

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.