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

: 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, ﬂowrates, 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 ﬂows 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 ﬂow was the most sensitive among all tested parameters. According to the Nash–Sutcli ﬀ e e ﬃ ciency (NSE), the MM_19 model is more e ﬃ cient 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 e ﬃ cient 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.


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 Water 2020, 12 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.
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 −1 ), 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.

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. where: -Mi = mass of PAH in compartment i; -Zi = fugacity capacity of compartment i; -fi and fj = fugacity of PAH in compartments i and j (Pa); -Vi = volume of compartment i (m 3 ); -Dji = transport D values to compartment i from compartment j with fugacity fj, including the flow in xylem and phloem and the uptake from soil and air (mol·Pa −1 ·h −1 ); -Do = transport and transformation D values for processes by which pollutant is removed from a giving compartment (mol·Pa −1 ·h −1 ).
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]. The process rate, N (mol·h −1 ), 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.

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. where: f i and f j = fugacity of PAH in compartments i and j (Pa); -V i = volume of compartment i (m 3 ); -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 −1 ·h −1 ); -D o = transport and transformation D values for processes by which pollutant is removed from a giving compartment (mol·Pa −1 ·h −1 ).
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].

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: and the mass balance of the roots compartment as follows: where: Replacing each D value and dividing Equations (2) and (3) by V L Z L and V R Z R , respectively, gives the following system of 2 coupled linear differential equations. 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: Analytical resolution of Equation (5) is provided in Supplementary Material A (SM_A).
Replacing each expression of C 1 and C 2 (SM_A), the analytical expressions of f L and f R are: All new constants are defined in SM_A.

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 −3 ) through a fugacity capacity Z (mol·m −3 ·Pa −1 ), such as: 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: PAH concentrations in stem_leaves (C L ) and roots (C R ) compartments can then be calculated as Z L f L and Z R f R , respectively. Rearranging Equation (9) and relating all fugacity capacity coefficients to water (W) gives 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.

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.

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 ( 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).
With %CO = %MO/1.72 and K OC = 0.411 K OW (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: where: Referring to Equation (8), K EW = 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: 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 −2 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 τ ij = ln (2) 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.

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: 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 −3 h −1 . The metabolic rate constant (k M ) was also assumed to be 2 × 10 −3 h −1 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  [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 −6 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.
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 2 ); (ii) a relatively humid period (average humidity = 70 ± 3%) and rainy (20.22 mm) with the highest average sunshine (220 ± 42 Wat/m 2 ); (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.

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 0 , 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.
Water 2020, 12, 1759 9 of 16 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.
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.

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.

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 −1 ) 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.

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.

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).

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.
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 KOW > 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)

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) showed that the model outputs were very sensitive to phloem flow and therefore to the multiplying factor (α). Table 1    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) 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 −5 ) for all species compared to the calibrated values of the Mackay_97 model.  Table 2 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 −1 -8.60 × 10 −1 ) in four cases (root compartment of E. indica and C. dactylon; stem_leaves compartment of A. sessilis and C. dactylon); satisfactory (5.50 × 10 −1 ) 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 −1 ); and very good (7.87 × 10 −1 ) just to simulate the soil-plant transfer of PAHs inside the root compartment of A. sessilis (Table 2).

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.
PAHs. Although the alpha factor of all plant species showed the same orders of magnitude, that of E. indica (6.64 × 10 −5 ) is the most reliable, since it was used for the other two species, with no significant differences between the initial and calibrated values.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/12/6/1759/s1, Table S1: Flowrates of PAHs from one compartment to another, Table S2: PAHs retention times in plant compartments, Table S3: PAHs half-lives of model compartments, Table S4: Example of sensitivity analysis for E. indica, Table S5: Sum of the square deviations (SSD) calculated between the simulated and the measured values for the 13 PAHs.