Modelling PAHs Transfer from Polluted Soil to Herbaceous Species in Phytoremediation Attempts

Abstract

To address the soil–plant transfer modelling of 13 US-EPA Polycyclic Aromatic Hydrocarbons (PAHs), a mechanistic model—MM_19—has been developed based on the fugacity concept. For that, the Mackay_97 model has been improved in terms of reconsidering the losses related to the transport and transformation mechanisms taking place in the compartments—roots and aboveground shoots—of the three short-life species (Eleusine indica, Cynodon dactylon and Alternanthera sessilis). Model input parameters consist of both experimental and literature data, including the initial soil and air PAHs content, flowrates, PAHs physico-chemical properties, retention times and transport half-lives of PAHs inside plant species. Using in situ weather data and Penman’s law, xylem flows were estimated as the evapotranspiration for each plant. Model calibration was performed using a Generalized Reduced Gradient (GRG) nonlinear optimization solver method. Sensitivity analysis showed that the phloem flow was the most sensitive among all tested parameters. According to the Nash–Sutcliffe efficiency (NSE), the MM_19 model is more efficient than the Mackay_97 model for all three plant species. Finally, the impact of PAHs physico-chemical parameters on their sol-plant transfer was discussed in terms of slight, intermediate and high molecules weight. The NSE values showed that the MM_19 model is more efficient than the Mackay_97 model. Indeed, comparisons between experimental and simulated results in the MM_19 model showed similarities for each compartment of the plant species. Thus, the MM_19 model can be used to predict the soil–plant transfer of organic pollutants.

1. Introduction

Soil–plant transfers of organic pollutants depend not only on the physico-chemical properties of the molecules under consideration but also on the physiological processes involved in plant development [1]. The goal of modelling organic pollutants in plants is to predict the kinetics and balance (of the extraction of these pollutants by species at different times) [2,3]. Both empirical and mechanistic models are usually used. Empirical models are based on statistical relationships and do not consider the physiological characteristics of the plant, nor the exchanges between compartments of the environment [4]. Thus, bioconcentration factors, regression equations, root concentration factor, partition coefficients and translocation factors have been established to predict the transfer of organic pollutants from soil to plant compartments [3]. There are some dynamic mechanistic models which integrate the physiological characteristics of plants, as well as the exchanges that can occur between plant and their environment (air, soil and/or water). There are different types of models, such as the one developed by [5] to predict the relationships between the lipophilicity, root uptake and translocation of non-ionised chemicals by barley. The first fugacity model that included a plant compartment was made by [6]. Then, the Trapp and Mackay group developed several models based on the same fugacity concept, including the Mackay_97 model based on three mass balance equations [2]. These models are all based on the same principles and describe the same basic processes such as advective absorption in plants, diffuse absorption, chemical equilibrium, xylem and phloem transport, growth dilution and deposition of soil and air particles [3]. Several studies addressed the soil–plant transfer modelling of organic pollutants [7] including PAHs [8]. These authors constructed an empirical model using two linear regression equations to express the concentration of 16 PAHs in lettuce, potatoes and carrots as a function of their concentrations in soil, without taking into account the possible impact of other soil characteristics.

To model the soil–plant transfer of PAHs, the Mackay_97 model is chosen in this paper because: (i) it is based on a compartmental approach that considers the dynamics of PAHs transfer through soils, roots, stems, leaves and air, (ii) the resulting differential equation could be solved analytically avoiding errors from numerical methods, (iii) this model applied to herbaceous plants with short vegetative periods, considers the physico-chemical properties of soils, the plant physiology and the key transfer processes between the plants compartments. The main objective of this work is to develop an improved version of the Mackay_97 model and use it to predict the transfer of 13 PAHs into 2 compartments (aerial and root) of 3 herbaceous plants as a function of time and initial PAH concentrations in soil and air. This will mainly involve: (i) improving the Mackay_97 model; (ii) conducting a sensitivity analysis and calibration of the model with the most sensitive parameter; (iii) assessing the efficiency of the two models; and (iv) comparing the simulated values with the measured ones.

2. Materials and Methods

2.1. Conceptual Model

The mechanistic model we are developing in this study (MM_19 model) takes into account 13 of the 16 US-EPA PAHs. The latter are those mainly analyzed inside the plants during the experimental study conducted by [9]. The 13 PAHs are fluorene (F), phenanthrene (Phe), anthracene (Ant), fluoranthene (Flu), pyrene (Pyr), benzo(a)anthracene (BaA), chrysene (Chy), benzo(b)fluoranthene (BbF), benzo(k)fluoranthene (BkF), benzo(a)pyrene (BaP), dibenzo(ah)anthracene (DahA), benzo(ghi)perylene (BP) et indeno(1,2,3-cd) pyrene (IP).

Two compartments are considered for each plant (Figure 1): aboveground shoots (L) = stems + leaves and roots (R). The other compartments are soil (S) and air (A). The main transfer flows across the plant are xylem (X) (upflow) and phloem (P) (downflow). Two keys processes within the compartments of the plant are taken into account: the growth dilution (G) and the metabolism (M) following the various degradation (catabolism) and synthesis (anabolism) reactions taking place inside the plant (Figure 1). According to [2], the Mackay_97 model considers the mass balance of downward and upward flows separately, hence avoiding a coupled equation. This leads to a poor consideration of losses related to the transport and transformation mechanisms that take place in each plant compartment; such as growth dilution and metabolism of PAHs inside the plants. Although approximate, this approach has the advantage to simplify the mathematical resolution of the Mackay_97 model. In order to overcome this, the present study is based on the development of a model that considers the mass balance of sap flow via the xylem and phloem vessels to be continuous, so that a single differential equation taking into account the two flows together could be established. Consequently, the mathematical resolution has been adapted. Another reason of improving the Mackay_97 model in this paper is due to the fact that it has not been calibrated.

The process rate, N (mol·h⁻¹), in or from a phase can be calculated by the product of transport (D) and fugacity (f). The straight arrows represent the process rate between compartments and the deformed arrows refer to the process rate of losses related to growth dilution and metabolism.

2.2. Model Equations

According to [2], for each compartment i of the plant, Equation (1) expressed the mass balance equation that governs the dynamics of pollutants across t compartments.

$$\frac{d\left(M_i\right)}{dt}=\frac{d\left(V_i{Z}_i{f}_i\right)}{dt}={{\displaystyle \sum }}^{ }{D}_{ji}{f}_j-f_i{{\displaystyle \sum }}^{ }{D}_o$$

(1)

where:

-

M_i = mass of PAH in compartment i;

-

Z_i = fugacity capacity of compartment i;

-

f_i and f_j = fugacity of PAH in compartments i and j (Pa);

-

V_i = volume of compartment i (m³);

-

D_ji = transport D values to compartment i from compartment j with fugacity f_j, including the flow in xylem and phloem and the uptake from soil and air (mol·Pa⁻¹·h⁻¹);

-

D_o = transport and transformation D values for processes by which pollutant is removed from a giving compartment (mol·Pa⁻¹·h⁻¹).

Transport and transformation processes are expressed in terms of D values that can be considered as a fugacity rate constant or a type of conductivity or the reciprocal of resistance. They are obtained from quantities such as flowrates, mass transfer coefficients, diffusivities and reaction rate constants [2].

2.2.1. Fugacity Expression of the Plant Compartments

Fugacity can be imagined as a measure of the escaping tendency of molecules from one phase to an adjacent phase. Hence, in a multicomponent system, if the fugacity of a component in the two adjacent phases is the same, the two phases will be in equilibrium with no net transfer of molecules from one phase to another. The fugacity can also be evaluated by simple methods, employing limiting assumptions, but still useful for engineering purposes at a variety of conditions [10]. The fluid and gaseous exchanges within the plant are governed by two flows (upward and downward).

According to Equation (1), PAH uptake terms (∑ D_ji f_j) include both the upwards and downwards input processes; e.g., for the root compartment, PAH uptake includes the stem_leaves to roots (D_LR) and the soil to roots (D_SR) transport processes. PAH output terms (f ∑D_o) also include both processes, e.g., the output terms for root compartment in both processes include roots to stem_leaves (D_RL) and root to soil transport (D_RS), along with metabolism (D_RM) and growth dilution (D_RG).

The resolution of Equation (1) leads to the expressed PAH fugacity and deduces its concentration inside the compartment under consideration according to the two flows direction. Thus, the mass balance of the stem_leaves compartment (L) can be written as:

$$\begin{gathered} \frac{d\left(V_L{Z}_L{f}_L\right)}{dt}=D_{AL}{f}_A+D_{RL}{f}_R-f_L{{\displaystyle \sum }}^{ }{D}_o \\ \mathrm{With} {{\displaystyle \sum }}^{ }{D}_o= D_{LA}+D_{LR}+D_{LM}+D_{LG} \end{gathered}$$

(2)

and the mass balance of the roots compartment as follows:

$$\begin{gathered} \frac{d\left(V_R{Z}_R{f}_R\right)}{dt}=D_{LR}{f}_L+D_{SR}{f}_S-f_R{{\displaystyle \sum }}^{ }{D}_o \\ \mathrm{With} {{\displaystyle \sum }}^{ }{D}_o=D_{RL}+D_{RS}+D_{LM}+D_{LG} \end{gathered}$$

(3)

where:

-

A = Air, R = Roots, S = Soil and L = Stems and Leaves

-

f_A, f_R, f_S and f_L are the fugacity of PAH in compartment A, R, S and L respectively;

-

D_AL, D_RL, D_SR and D_LR are the air to stem_leaves, roots to stem_leaves, soil to roots and stem_leaves to roots transport D values of PAH. With $D=\ln \left(2\right)\left(VZ\right)/\tau$, with $\tau$ the PAH half-life.

-

D_LM, D_LG, are the PAH lost by metabolism and growth dilution. D_LA, D_LR,D_RS and D_RL are the other losses related to the first compartment. For example, D_LA is the transport D value lost from stem_leaves to air (through transpiration for instance).

Replacing each D value and dividing Equations (2) and (3) by V_LZ_L and V_RZ_R, respectively, gives the following system of 2 coupled linear differential equations.

$$\{\begin{array}{l}\frac{d{f}_L}{dt} = -\left[\frac{\ln \left(2\right)}{{\tau }_{LA,L}}+\frac{\ln \left(2\right)}{{\tau }_{LR,L}} + k_{LG} + k_{LM}\right]f_L + \frac{\ln \left(2\right)}{{\tau }_{RL,L}}{f}_R + \frac{\ln \left(2\right)}{{\tau }_{AL,L}}{f}_A\\ \frac{d{f}_R}{dt} = -\left[\frac{\ln \left(2\right)}{{\tau }_{RL,R}} + \frac{\ln \left(2\right)}{{\tau }_{RS,R}} + k_{RG} + k_{RM}\right]f_R + \frac{\ln \left(2\right)}{{\tau }_{LR,R}}{f}_L + \frac{\ln \left(2\right)}{{\tau }_{SR,R}}{f}_S\end{array}$$

(4)

where:

-

τ_AL,L, τ_RL,L, τ_LA,L and τ_LR,L are the air to stem_leaves, root to stem_leaves, stem_leaves to air and stem_leaves to root PAH clearance half-lives in the stem_leaves compartment;

-

τ_LR,R, τ_SR,R, τ_RL,R and τ_RS,R are the stem_leaves to root, soil to root, root to stem_leaves and root to soil PAH clearance half-lives in the root compartment;

-

k_LG and k_LM are stem_leaves growth dilution and metabolism rates constants, respectively; k_RG and k_RM the roots growth dilution and metabolism rates constants, respectively.

Equation (4) can be written as a linear differential equation with the second element:

$$\begin{gathered} \frac{df}{dt}=Mf+F \\ \mathrm{where}: f=\left(\begin{array}{c}{f}_L\\ f_R\end{array}\right); F=\left(\begin{array}{c}\frac{\ln \left(2\right)}{{\tau }_{AL,L}}{f}_A\\ \frac{\ln \left(2\right)}{{\tau }_{SR,R}}{f}_S\end{array}\right) \mathrm{and} \\ M=\left[\begin{array}{cc}-\left[\frac{\ln \left(2\right)}{{\tau }_{LA,L}}+\frac{\ln \left(2\right)}{{\tau }_{LR,L}}+k_{LG}+k_{LM}\right]& \frac{\ln \left(2\right)}{{\tau }_{RL,L}}\\ \frac{\ln \left(2\right)}{{\tau }_{LR,R}}& -\left[\frac{\ln \left(2\right)}{{\tau }_{RL,R}}+\frac{\ln \left(2\right)}{{\tau }_{RS,R}}+k_{RG}+k_{RM}\right]\end{array}\right] \end{gathered}$$

(5)

Analytical resolution of Equation (5) is provided in Supplementary Material A (SM_A).

Replacing each expression of C₁ and C₂ (SM_A), the analytical expressions of f_L and f_R are:

$$\begin{array}{ll}{f}_L & = \left[\frac{\left(A+{\lambda }_2\right)}{\left({\lambda }_1-{\lambda }_2\right)}\left[\frac{\ln \left(2\right)}{{\tau }_{AL,L}}{f}_A - \frac{\mathrm{Dln}\left(2\right)}{\left(A+{\lambda }_2\right){\tau }_{SR,R}}{f}_S\right]\left(\frac{1}{{\lambda }_1}\right)\left(1 - e^{{\lambda }_{1}t}\right)\right]\\ & + \left[\frac{\left(A+{\lambda }_1\right)}{\left({\lambda }_2-{\lambda }_1\right)}\left[\frac{\ln \left(2\right)}{{\tau }_{AL,L}}{f}_A - \frac{\mathrm{Dln}\left(2\right)}{\left(A+{\lambda }_1\right){\tau }_{SR,R}}{f}_S\right]\left(\frac{1}{{\lambda }_2}\right)\left(1 - e^{{\lambda }_{2}t}\right)\right]\end{array}$$

(6)

$$\begin{array}{ll}{f}_R =& \left[\frac{\left(A+{\lambda }_1\right)}{D} \times \frac{\left(A+{\lambda }_2\right)}{\left({\lambda }_1-{\lambda }_2\right)}\left[\frac{\ln \left(2\right)}{{\tau }_{AL,L}}{f}_A - \frac{\mathrm{Dln}\left(2\right)}{\left(A+{\lambda }_2\right){\tau }_{SR,R}}{f}_S\right]\left(\frac{1}{{\lambda }_1}\right)\left(1 - e^{{\lambda }_{1}t}\right)\right]\\ + \left[\frac{\left(A+{\lambda }_2\right)}{D} \times \frac{\left(A+{\lambda }_1\right)}{\left({\lambda }_2-{\lambda }_1\right)}\left[\frac{\ln \left(2\right)}{{\tau }_{AL,L}}{f}_A - \frac{\mathrm{Dln}\left(2\right)}{\left(A+{\lambda }_1\right){\tau }_{SR,R}}{f}_S\right]\left(\frac{1}{{\lambda }_2}\right)\left(1 - e^{{\lambda }_{2}t}\right)\right]\end{array}$$

(7)

All new constants are defined in SM_A.

2.2.2. PAHs Concentrations in the Roots and Stem_Leaves Compartments

Fugacity (f) which is an equilibrium criterion and can be considered as partial pressure (Pa), is linearly related to concentration C (mg·kg⁻³) through a fugacity capacity Z (mol·m⁻³·Pa⁻¹), such as:

$$C=fZ$$

(8)

According to [2], equilibrium partitioning (at equal fugacity, f_i = f_j) between compartment, with concentrations C_L and C_R, can be described by a dimensionless partition coefficient (K_ij) which is essentially a ratio of Z values; hence:

$$K_{ij}=\frac{C_i}{C_j}=\frac{f_i{Z}_i}{f_j{Z}_j}=\frac{Z_i}{Z_j}$$

(9)

PAH concentrations in stem_leaves (C_L) and roots (C_R) compartments can then be calculated as Z_Lf_L and Z_Rf_R, respectively. Rearranging Equation (9) and relating all fugacity capacity coefficients to water (W) gives

$$\begin{array}{ll}{C}_L =& \left[\frac{\left(A + {\lambda }_2\right)}{\left({\lambda }_1 - {\lambda }_2\right)}\left[\frac{\ln \left(2\right)}{{\tau }_{AL,L}}{C}_A\left(\frac{K_{LW}}{K_{AW}}\right)- \frac{\mathrm{Dln}\left(2\right)}{\left(A+{\lambda }_2\right){\tau }_{SR,R}}{C}_S\left(\frac{K_{LW}}{K_{SW}}\right)\right]\left(\frac{1}{{\lambda }_1}\right)\left(1 - e^{{\lambda }_{1}t}\right)\right]\\ + \left[\frac{\left(A+{\lambda }_1\right)}{\left({\lambda }_2-{\lambda }_1\right)}\left[\frac{\ln \left(2\right)}{{\tau }_{AL,L}}{C}_A\left(\frac{K_{LW}}{K_{AW}}\right) -\frac{\mathrm{Dln}\left(2\right)}{\left(A+{\lambda }_1\right){\tau }_{SR,R}}{C}_S\left(\frac{K_{LW}}{K_{SW}}\right)\right]\left(\frac{1}{{\lambda }_2}\right)\left(1 - e^{{\lambda }_{2}t}\right)\right]\end{array}$$

(10)

$$\begin{array}{ll}{C}_R& = \left[\frac{\left(A + {\lambda }_1\right)}{D}\times \frac{\left(A + {\lambda }_2\right)}{\left({\lambda }_1 - {\lambda }_2\right)}\left[\frac{\ln \left(2\right)}{{\tau }_{AL,L}}{C}_A\left(\frac{K_{RW}}{K_{AW}}\right)- \frac{\mathrm{Dln}\left(2\right)}{\left(A + {\lambda }_2\right){\tau }_{SR,R}}{C}_S\left(\frac{K_{RW}}{K_{SW}}\right)\right]\left(\frac{1}{{\lambda }_1}\right)\left(1 - e^{{\lambda }_{1}t}\right)\right]\\ & + \left[\frac{\left(A + {\lambda }_2\right)}{D} \times \frac{\left(A + {\lambda }_1\right)}{\left({\lambda }_2 - {\lambda }_1\right)}\left[\frac{\ln \left(2\right)}{{\tau }_{AL,L}}{C}_A\left(\frac{K_{RW}}{K_{AW}}\right)-\frac{\mathrm{Dln}\left(2\right)}{\left(A + {\lambda }_1\right){\tau }_{SR,R}}{C}_S\left(\frac{K_{RW}}{K_{SW}}\right)\right]\left(\frac{1}{{\lambda }_2}\right)\left(1 - e^{{\lambda }_{2}t}\right)\right]\end{array}$$

(11)

where K_LW, K_AW, K_SW and K_RW are the stem_leaves to water, air to water, soil to water and root to water partition coefficients, respectively.

2.3. Model Input Parameters

The input parameters of the model were either measured or calculated using experimental or literature data, or were taken directly from the literature without any modifications.

The calculated input parameters are:

-

Partition coefficients of PAHs between compartments;

-

Flowrates that control the transport processes;

-

Retention times of PAHs in each compartment of the plant according to the flow direction;

-

Half-lives of transport (τ), growth dilution (τ_G) and metabolism (τ_M).

The existing input parameters are:

-

Metabolism and growth dilution half-lives

-

Plant morphological and physiological parameters such as volumes of stem_leaves and roots compartments for each plant species; and the initial soil PAHs concentrations.

2.3.1. Calculated Input Parameters

Partition Coefficients and the Initial PAHs Concentration in Water

To calculate the initial PAHs concentration in water (C_W), the distribution coefficient (K_d = C_S/C_W) is used. Thus:

$$C_W=C_S/K_d,$$

(12)

where C_S is the initial soil PAH concentration (mg/kg).

K_d is calculated using the organic carbon coefficient (K_OC: ratio between the amount of PAH absorbed per unit weight of organic carbon and the amount of PAH in the aqueous phase at equilibrium).

$$\begin{gathered} K_d=\frac{K_{OC} \times \%CO }{100} \\ \mathrm{With} \%CO=\%MO/1.72 \mathrm{and} K_{OC}=0.411 K_{OW} \end{gathered}$$

(13)

where %CO [11] and %MO are the percentages of organic carbon and organic matter, respectively, from experimental data.

The compartment-water partition coefficient (K_iW) including stem_leaves, root and soil, have been calculated using the leaves (E) to air (A) partition coefficient (K_EA) from [12] as follows:

$$\begin{gathered} K_{EA}=\left(\raisebox{1ex}{$RT$}\!\left/ \!\raisebox{-1ex}{$H$}\right.\right)\times \left[\begin{array}{c}{V}_{C}1.11K_{OW}^{0.97}+V_W+V_l{K}_{OW}+V_{F}0.0372K_{OW}^{0.95}+V_P(86.2K_{OW}^{-1}\\ +3.70)\end{array}\right] \\ \mathrm{With}: \raisebox{1ex}{$RT$}\!\left/ \!\raisebox{-1ex}{$H$}\right.=K_{AW} \end{gathered}$$

(14)

where:

-

V_C, V_W, V_l, V_F and V_P are the volume fractions of cuticular membrane, water, carbohydrate, lipids and protein of the leaves.

-

K_AW is the air to water partition coefficient and RT/H derived from the Henry’s law constant.

Referring to Equation (8), $K_{EW}=\frac{K_{EA}}{K_{AW}}=K\_iW$. Rearranging Equation (14) and taking into account the density of each compartment i, the final expression of K_iW is then:

$$K_{iW}=\left[V_{C}1.11K_{OW}^{0.97}+V_W+V_L{K}_{OW}+V_{F}0.0372K_{OW}^{0.95}+V_P\left(86.2K_{OW}^{-1}+3.70\right)\right]\times d_i$$

(15)

The octanol/cuticle-air partition coefficient can be expressed as K_OA = K_OW/K_AW.

Flowrates and Retention Times

Regarding the case of this study with two plant compartments (stem_leaves = L and root = R), the calculation of flowrates related to stem_leaves and air compartments (Table S1) is deduced from [2] and all the explanations related to the transport parameters are in Supplementary Material B (SM_B). Flowrates are calculated according to the migration direction of PAHs. Xylem flowrates (G_X) are specific to the plant type and dependent on weather conditions. According to [2], the Phloem flow is equal 5% of the xylem flow (G_P = 5 × 10⁻² G_X); the soil to root flowrate is 1.05 of the xylem flow (G_SR = 1.05 G_X); and the root to soil flowrate is 5% (α) of the xylem flow (G_RS = 0.05 G_X).

Retention times (𝜹)—the times required for a volume of xylem or phloem sap to equal the volume of a given compartment to flow through the plant—are calculated from flowrates as presented in Table S2 (SM_B).

Transport Half-Lives

Half-lives D values for each flow direction are expressed in Table S3 (Supplementary Material C). It is the time required to halve the amount of PAH in a given compartment. The transport half-lives (D = ln(2)VZ/τ) are deduced from the two plant compartments volumes (V_L and V_R) and the four flowrates (G_LA, G_AL, G_P and G_SR). According to [2], the general expression to calculate the clearance half-life of a PAH from a compartment i to a compartment j, G_ij is ${\tau }_{ij}=\ln \left(2\right)K_{iw}\frac{V_i}{G_{ij}}=\ln \left(2\right)K_{iw}{\delta }_{i,ij}.$ Thus, half-lives are then calculated for each flow direction using the values of retention times (𝜹) and the partition coefficients.

To determine the air–leaves exchanges, the transfer half-lives are explained in Supplementary Material C.

2.3.2. Existing Input Parameters

Metabolism and Growth Dilution Half-Lives

Most persistent pollutant of environmental concern like PAHs have low metabolic rates and the metabolism half-life can be assumed to be infinite or very large [2]. The D values for metabolism and growth dilution, for example, in the stem_leaves compartment, are D_LM and D_LG, respectively. These D values can be calculated as:

$$D_{LG}=Z_L\left(\frac{d{V}_L}{dt}\right)=k_{LG}{Z}_L{V}_L$$

(16)

$$D_{LM}=Z_L\left(\frac{d{V}_L}{dt}\right)=k_{LM}{Z}_L{V}_L$$

(17)

where V is the volume of the plant, and k_LM and k_LG are the assumed first order rate constants. τ_LM and τ_LG, the stem_leaves compartment metabolic and growth dilution half-life are ln(2)/k_LM and ln(2)/k_LG respectively.

Later, as the plant matures, k_LG decreases. For simplicity, it is assumed that k_LG is constant throughout the short growing period of an herbaceous agricultural plant. The growth half-life, τ_LG, is usually the simple estimation of the time required for each plant tissue to double in volume. For grasses, τ_LG may be a week [2].

According to experimental data, the time required for the three plant species to double their original size from sprouting is approximately 14 days = 336 h [9]; corresponding to a k_G of 2 × 10⁻³ h⁻¹. The metabolic rate constant (k_M) was also assumed to be 2 × 10⁻³ h⁻¹ for all plant compartments, corresponding to a metabolic half-life of 336 h.

Plant Morphological and Physiological Parameters and Others

Mass, volumes and density of each compartment (air, stem_leaves, roots, soil and water) were obtained from the experimental data and literature. The expression of PAH content in the atmosphere is developed in Supplementary Material D (SM_D). Volume fractions of cuticular membrane, carbohydrate, lipids and protein of each plant compartment come from [12], and the volume fraction of water in each plant compartment has been calculated using their water content. The typical values for τ_O and τ_A are 126 h and 5 × 10⁻⁶ h, respectively, and the diffusion fraction of xylem flow (Ø) is 0.05 [2]. The duration used to model the transfer of PAHs (150 days = 3750 h) is the same as that of the experiments conducted by [9]. The octanol-water partition coefficient (log K_OW) of each PAH comes from literature as well as their Henry law values. All these input parameters were used to solve the model.

2.4. Overview of the Experimental Work

An experimental study has been conducted by [13] to investigate the effect of fuel oil pollution on growth characteristics and phytoremediation potential of six plant species (Cynodon dactylon (L.) Pers., Eleusine indica (L.) Gaertn., Cleome ciliata Schum. and Thonn., Alternanthera sessilis (L.) R. Br. ex DC†, Commelina benghalensis L. and Asystasia gangetica (L.) T. Anderson).

The experimental device comprises three modalities, namely Tn: unpolluted planted soil; To: unplanted polluted soil; and Tp: polluted planted soil. For this study, soil had sandy loam texture (81.06% sand, 11.25% silt, and 3.98% clay); with a low cation exchange capacity (10.38 ± 1.08 meq/100 g of soil), a pH (7.4 ± 0.26) close to the neutrality; rich in organic matter (5.86 ± 0.51% DM) and nutrients (total Kjeldahl nitrogen = 1516.67 ± 44.1 mg/kg DM; total phosphorus = 1326.51 ± 31.79 mg/kg DM). The study was carried out under natural conditions at an experimental site located in the North-East of Strasbourg city (48°34′24.21′ N; 7°45′8.47′ E), France. A WatchDog weather station was used for climate data registration. The climatic variations were diverse because the experimental period (April–August) took place over two climatic seasons (spring and summer). A hierarchical classification makes it possible to distinguish three climates: (i) a fairly cold period (average temperature = 14 ± 1 °C) and wet (average humidity = 67 ± 5%), with the lowest average sunshine (153 ± 31 Wat/m²); (ii) a relatively humid period (average humidity = 70 ± 3%) and rainy (20.22 mm) with the highest average sunshine (220 ± 42 Wat/m²); (iii) and a warm period (average temperature = 24 ± 1 °C) and dry (average humidity 59 ± 1%).

As for the results of this experiment, we noticed that only 3 plants (E. indica, C. dactylon and A. sessilis) survived and developed throughout the 150 days of experimentation. The bioconcentration and translocation factors values of PAH indicate that E. indica and A. sessilis promoted the rhizodegradation and phytoextraction of hydrocarbon-polluted soils, whereas C. dactylon was only involved into rhizodegradation. Input parameters from this study are the initial PAHs concentration in soil; the mass and density of compartment R, L and S; and the weather data to calculate plants transpiration (SM_B) used as xylem flow from which the phloem flow was deduced.

2.5. Model Resolution and Data Analysis

The two linear differential equations were solved for the 13 PAHs and each of the 3 plant species (Eleusine indica, Cynodon dactylon and Alternanthera sessilis).

Sensitivity analysis was mainly based on the phloem flow calculation factor (α), the diffusion fraction of xylem flow between plant compartments (Ø), as well as the characteristic time of PAH through cuticle (𝝉_O). The choice of these parameters is based on the fact that they come from literature and are likely to vary according to the type of pollutant and the environmental factors. They are also directly related to the plant.

The classical Mean Square Error (MSE) here measures the difference between computed values h and the measured values h⁰, for n PAHs. MSE is the one of most widely used criteria for calibration and evaluation of models with observed data. It can be expressed in the units of the quantities it compares when a unique variable is considered.

$$\mathrm{MSE}=\frac{1}{n} {\displaystyle \sum }_{i=1}^n{\left({\mathrm{h}}_{\mathrm{i}}-{ \mathrm{h}}_{\mathrm{i}}^0\right)}^2 \in \left(0,+\infty \right)$$

(18)

The Generalized Reduced Gradient (GRG) nonlinear optimization solver method [14] in Excel was used for the calibration.

To assess model efficiency, the so-called Nash–Sutcliffe efficiency (NSE), a dimensionless error criterion, has been used. It accounts for the mean μ^o and standard deviation σ^o of the dataset.

$$\mathrm{NSE}=1-\frac{ {{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{n}}{\left({\mathrm{h}}_{\mathrm{i}}-{\mathrm{h}}_{\mathrm{i}}^0\right)}^2 }{{{\displaystyle \sum }}_{\mathrm{t}=1}^{\mathrm{n}}{\left({\mathrm{h}}_{\mathrm{i}}^0-{\mu }^0\right)}^2 }=1-\frac{\mathrm{MSE}}{{\left({\sigma }^0\right)}^2}\in \left(-\infty , 1\right)$$

(19)

A value of 1 for NSE indicates a perfect agreement between the model and the dataset. Model performance is rated as “very good” if 0.75 < NSE < 1.00, as “good” if 0.65 < NSE < 0.75, as “satisfactory” if 0.50 < NSE < 0.65, and “unsatisfactory” otherwise.

The Student’s t-test was used to compare the average values of PAHs contents between the computed and the measured values; using StatPlus:mac LE software v6.1.60 (AnalystSoft Inc., Walnut, California). Evapotranspiration (ETP) was calculated using the Penman–Monteith method.

3. Results and Discussion

3.1. Sensitivity Analysis

The sensitivity analysis shows that changes in the α factor have the greatest influence on the model (Table S4, Supplementary Material E). In fact, changing the α factor by an order of magnitude results in equivalent changes in the simulated values. Parameter 𝝉_O also significantly influences the model; however, its impact is not homogeneous like the α factor. As for Ø, the model is less sensitive to its variation compared to the other two parameters. Indeed, it only influences the root compartment.

For example, in the case of chrysene (Chy) in the stem_leaves compartment, by multiplying alpha by 10⁻¹ and 10⁻², the initial value of Chy (4.44 × 10⁻¹) is modified by 10⁻¹ (4.59 × 10⁻²) and 10⁻² (4.63 × 10⁻³), respectively. Applying the same scenario with 𝝉_O, no change is observed when 𝝉_O is multiplied by 10⁻¹ (4.47 × 10⁻¹); and when 𝝉_O is multiplied by 10⁻², the simulated Chy values is modified in the order of 10⁻¹ (3.67 × 10⁻¹). Multiplying Ø by 10⁻², the initial Chy value (7.06 × 10⁻³) in the roots is only modified by 10⁻¹ (6.75 × 10⁻³) (Supplementary Material E).

3.2. Impact of Light PAHs on Model Calibration

Model calibration with the 13 PAHs was not possible; it was performed with 10 PAHs. The latter are the heaviest PAHs (4–6 benzene rings). The remaining 3 PAHs are 3-ring (Fluorene, Anthracene and Phenanthrene), thus the lower molecular weight PAHs among the 13 US-EPA PAHs used in this study. Overall, whatever the plant compartment and the species, the sum of the square deviations (SSD) calculated between the simulated and the measured values for these 3 PAHs are not in the same orders of magnitude as those of the other 10 PAHs (Table S5, Supplementary Material F). Thus, this contributes to skewing the model calibration.

Different behaviour observed between heavy and light PAHs during modelling could be explained by their physico-chemical properties as well as the environment factors. Regarding PAHs properties, it is assumed that soil–plant transfer of PAHs is only possible for low molecular weight PAHs (2 to 3 benzene rings) while those of 3 to 6 rings show a tendency towards adsorb to the cortical zone of roots [15]. In addition, light PAHs (152 to 178 g·mol⁻¹) have higher solubility and volatility than heavy PAHs which increase their probability to be transfer into the plant through xylem and phloem flows. It is at this stage that environmental factors are involved. In fact, for experiments, with regard to phloem flow, the physical non-restriction of the air compartment and thus the increase in volatilization would reduce downward transfers; and for xylem flows, the soil biodegradation of these highly soluble PAHs (2 to 3 benzene rings) must be taken into account, since plant–microbe interaction significantly affect the bioavailability and remediation of PAHs [16]. Moreover, the rhizosphere could also stabilize pollutants through polymerization reactions such as humification [17]. All these statements make it possible to understand why the model overestimated (Supplementary Material E) the transfer of 3-ring PAHs (Fluorene, Anthracene and Phenanthrene), especially within the root compartment, for almost all three herbaceous species. With respect to the foregoing, only 10 PAHs will be taken into account for both model calibration and modelling.

3.3. Impact of the Physico-Chemical Properties of PAHs on Their Transfer from Soil to the Plant Species

Overall, the soil–plant transfer modelling results of the heavier PAHs follow the same trend, regardless of the plant species. In fact, the simulated concentrations of heavy PAHs are usually closer to the experimental ones. However, some differences are noticeable depending on the physico-chemical properties of PAHs.

3.3.1. Intermediate Molecular Weights’ PAHs

Fluoranthene (Flu) and Pyrene (Pyr) are among the intermediate molecular weights’ PAHs [18]. Compared to the eight other PAHs, Flu and Pyr concentrations present some significant differences between the MM_19 model and the experimental results. In fact, the modelling results of Flu in the roots parts (Figure 2B,D,F) are significantly higher than the experimental results for all the three plant species (A. sessilis: p < 0.001, C. dactylon and E. indica: p < 0.05). In the stem_leaves compartment of C. dactylon, the simulated values of Pyr are significantly lower (p < 0.01) than their experimental ones (Figure 2F).

The higher Flu concentrations for the MM_19 model compared to the experiments predict a strong bioaccumulation trend for this PAH. Indeed, Flu is intermediate between the most soluble and volatile PAHs and the most bioaccumulative PAHs (log K_OW > 3). Thus, its tendency to be more or less soluble in the soil solution and to have a strong affinity with the plants lipids contributes to its bioaccumulation inside the plant roots [15].

The low simulated Pyr concentrations compared to the experimental ones can be explained by the adsorption phenomenon. The difference appears only for one plant species (C. dactylon), so it is assumed that the adsorbed fraction of Pyr on the roots of C. dactylon has not been desorbed, probably due to the high affinity that could exist between Pyr and the roots lipids of C. dactylon. In fact, PAHs are well-known to be lipophilic.

The Pyr and Flu contents of the Mackay_97 model are significantly different from those tested in some plant species. Indeed, in A. sessilis, the simulated values of Pyr and Flu in the stem_leaves compartment of the Mackay_97 model are significantly higher (p < 0.001) than those of the experiments (Figure 2B). The same observation is made for E. indica regarding Pyr concentration (p < 0.05) in the stem_leaves compartment (Figure 2A). No significant differences were observed for the root parts of each species (Figure 2B,D,F). The comparison between the simulated values of different models shows that the root concentrations of the Mackay_97 model are slightly different (p < 0.1) to the experimental ones compared to the MM_19 model (p > 0.05).

3.3.2. High Molecular Weights’ PAHs

Simulated (MM_19 model) and measured concentrations of the eight other heavy PAHs are close to the experimental concentrations, except the simulated root concentration of Benzo(a)anthracene (BaA) which is higher than the experimental one for A. sessilis (Figure 3D). Furthermore, the MM_19 model underestimates the experimental concentrations of Benzo(ghi)perylene (BP) (p < 0.001) and Indeno(1,2,3-cd)pyrene (IP) (p < 0.05) in the stem_leaves of all plant species (Figure 3A,C,E).

The close concentrations of heavy PAHs between the MM_19 model and the experiments in the two plant compartments of each herbaceous species confirm their low bioavailability as well as their high bioaccumulation following their log K_OW > 3. In fact, heavy PAHs tend to be accumulated on plants lipids and high organic matter area [19].

Comparing the Mackay_97 model with experiments shows that it underestimates the experimental concentrations of all the eight PAHs (p < 0.01) in the stem_leaves of A. sessilis (Figure 3C). The BP simulated concentration with the Mackay_97 model is the only one of the eight PAHs that is significantly higher than the experimental ones in the root parts of C. dactylon (p < 0.01) (Figure 3F) and E. indica (p < 0.05) (Figure 3B).

Overall, although the differences are not statistically perceptible, with 95% certainty between the experimental concentrations and those of the Mackay_97 model, the Student’s t-test shows differences with 90% certainty—particularly for Chy and BbF, whose experimental concentrations are higher than those of the Mackay_97 model for C. dactylon and E. indica. This is not the case with the MM_19 model. Thus, the experimental values of PAHs concentrations in the plant compartments are closer to the model developed in this study (the MM_19 model) compared to that of Mackay_97. Probably because of the difference between the two models in terms of the initial mass balance of the upward and downward flowrates.

3.4. Comparison Between the MM_19 Model and the Mackay_97 Model

3.4.1. Significance of the Phloem Flow Multiplying Factor (α)

As specified in the methodology part, the phloem is calculated from the xylem flow using a multiplying factor (α). The difficulty in estimating this multiplying factor is mainly due to the fact that the air compartment is not physically limited as the model considers it. Direct simulations carried out by varying the model’s input interest parameters (

Table 1. Initial and calibrated values of the phloem flow multiplying factor (α).

Models	Phloem Flow Multiplying Factor (α)	E. indica	C. dactylon	A. sessilis
MM_19	Initial value	5.00 × 10−2	6.64 × 10−5	6.64 × 10−5
	Calibrated value	6.64 × 10−5	6.95 × 10−5	9.22 × 10−5
Mackay_97	Initial value	5.00 × 10−2	2.72 × 10−4	2.72 × 10−4
	Calibrated value	2.72 × 10−4	6.07 × 10−4	4.69 × 10−2

) showed that the model outputs were very sensitive to phloem flow and therefore to the multiplying factor (α). Table 1 presents the values of α factor for all plant species, before and after calibration. The calibrated value obtained for E. indica was used as the initial value (before calibration) for C. dactylon and A. sessilis. There are significant differences between the phloem flow of the Mackay_97 model and that of the MM_19 model. In fact, bromacil, the pollutant used by [2], has low values of H and log Kow. It could therefore be easily transferred through phloem and xylem flows; hence the high factor of phloem flow. Moreover, it appears that the Mackay_97 model was not calibrated first. Table 1 also shows that the calibrated values of α factor used to estimate phloem flow for the MM_19 model have some consistency, given the same orders of magnitude (10⁻⁵) for all species compared to the calibrated values of the Mackay_97 model.

3.4.2. Models Efficiency

Table 2. Efficiency of MM_19 and Mackay_97 models.

NSE	Models	E. indica	C. dactylon	A. sessilis
NSE_L	MM_19	5.50 × 10−1	8.07 × 10−1	8.45 × 10−1
	Mackay_97	5.87 × 10−1	5.20 × 10−1	−6.56 × 103
NSE_R	MM_19	8.60 × 10−1	7.24 × 10−1	−1.05 × 102
	Mackay_97	−1.47 × 10−1	−1.95 × 10−1	7.87 × 10−1

provides the efficiency of the two models through the NSE values. Overall, it shows that the MM_19 model is more efficient than the Mackay_97 model. NSE values of MM_19 model show that its efficiency is very good (7.24 × 10⁻¹–8.60 × 10⁻¹) in four cases (root compartment of E. indica and C. dactylon; stem_leaves compartment of A. sessilis and C. dactylon); satisfactory (5.50 × 10⁻¹) for the stem_leaves compartment of E. indica; and unsatisfactory just for the root compartment of A. sessilis. The efficiency of Mackay_97 model is unsatisfactory in three cases (root compartment of E. indica and C. dactylon; stem_leaves compartment of A. sessilis); satisfactory for the stem_leaves compartment of E. indica (0.587) and C. dactylon (5.20 × 10⁻¹); and very good (7.87 × 10⁻¹) just to simulate the soil–plant transfer of PAHs inside the root compartment of A. sessilis (Table 2).

3.5. Modelling Results

The cluster dendrogram and the PCA (Figure 4) compare the simulated and experimental values of the 10 PAHs modelled for the MM_19 model. They do not take into account the PAH values in the root compartment of A. sessilis for which the effectiveness of the MM_19 model was unsatisfactory. These comparisons show in general the similarities between model and experiments for the two compartments of each plant. In fact, in each of the two cluster groups, the experimental and simulated values corresponding to a compartment of any plant are present. For example, in the first group, both experimental and simulated PAH values are found in the root compartment of E. indica, and C. dactylon.

Overall, the simulated concentrations of the 10 heavy PAHs are generally closer to the experimental concentrations for the different compartments of all plant species. However, some differences are noticeable depending on the physico-chemical properties of PAHs. More specifically, the model underestimated the experimental concentrations of four PAHs, mainly those of highest molecular weight (DahA, BP and IP) in the air compartments of all plants. This result is consistent with the literature since, due to their low bioavailability and high bioaccumulation, heavy PAHs tend to accumulate on plant lipids and areas with a high organic matter content [19]. Simulated Flu concentrations in the root parts are significantly higher than those in the experimental parts (p < 0.05) for C. dactylon and E. indica. This predicts a strong tendency for the bioaccumulation and transfer of this PAH to the roots. Only the simulated values of Pyr are significantly lower (p < 0.01) than the experimental values in the stem _leaves compartment of C. dactylon.

4. Conclusions

This work aimed to model the soil–plant transfer of 13 US-EPA PAHs using the fugacity concept of Mackay_97 model. However, the double integration of losses in this model, as well as its non-calibration, led to the development of a new one, the MM_19 model. Sensitivity analysis shows that the phloem flow multiplying factor (α) was the most sensitive parameter of the study. It has then been used for model calibration. Of the 13 PAHs, only the soil–plant transfer of 10 (4–6 benzene rings) was simulated due to significant differences in the sum of the square deviations found for the 3-ring PAHs (F, Ant and Phe), which caused a skewness in the model’s calibration. Analysing the performance of both models revealed that the MM_19 model is more efficient than the Mackay_97 model. Solving MM_19 model showed that it overestimated the measured values of two intermediate PAHs (Flu and Pyr); especially in the roots of all species for Flu, and in the stem_leaves of C. dactylon for Pyr. The other eight heaviest PAHs showed significant similarity between the model and measured values; except for DahA, BP, and IP, the measured values of which were underestimated by the model in the stem_leaves compartment of each species. According to its high efficiency for all the three herbaceous species and its accuracy in predicting the soil–plant transfer of the 8 on 10 PAHs, the MM_19 model can be recommended for similar work. However, the integration of rhizospheric biodegradation processes of PAHs by microorganisms would be useful to provide more accurate predictions of the MM_19 model. Results obtained from the calibrated values of the α factor in the MM_19 model make it possible to recommend its future use in the modelling process for the specific case of lipophilic pollutants like PAHs. Although the alpha factor of all plant species showed the same orders of magnitude, that of E. indica (6.64 × 10⁻⁵) is the most reliable, since it was used for the other two species, with no significant differences between the initial and calibrated values.