On the Use of Mechanistic Soil-Plant Uptake Models: A Comprehensive Experimental and Numerical Analysis on the Translocation of Carbamazepine in Green Pea Plants.

Food contamination is a major worldwide risk for human health. Dynamic plant uptake of pollutants from contaminated environments is the preferred pathway into the human and animal food chain. Mechanistic models represent a fundamental tool for risk assessment and the development of mitigation strategies. However, difficulty in obtaining comprehensive observations in the soil-plant continuum hinders their calibration, undermining their generalizability and raising doubts about their widespread applicability. To address these issues, a Bayesian probabilistic framework is used, for the first time, to calibrate and assess the predictive uncertainty of a mechanistic soil-plant model against comprehensive observations from an experiment on the translocation of carbamazepine in green pea plants. Results demonstrate that the model can reproduce the dynamics of water flow and solute reactive transport in the soil-plant domain accurately and with limited uncertainty. The role of different physicochemical processes in bioaccumulation of carbamazepine in fruits is investigated through Global Sensitivity Analysis, which shows how soil hydraulic properties and soil solute sorption regulate transpiration streams and bioavailability of carbamazepine. Overall, the analysis demonstrates the usefulness of mechanistic models and proposes a comprehensive numerical framework for their assessment and use.

Introduction

Food contamination from polluted environments represents a major risk for human health. Plants represent the most common pathway into the human and animal food chain for environmental pollutants. Chemicals are taken up by plant roots and then translocated toward the edible parts, where they are bioaccumulated[1−6] or eventually metabolized in byproducts,[7,8] which can be toxic.[9,10] Soil is among the most important pollution sources for plants.[11] It has been shown that it can be contaminated by human and veterinary pharmaceuticals,[12−14] uptake of which is controlled by plants, chemicals, and soil properties.[15−17]

Carbamazepine (CBZ) is a widely prescribed anticonvulsant pharmaceutical and antiepileptic drug. Once administered, it is degraded by the liver in multiple metabolites, which are then excreted in both feces and urine. Due to its persistence, it is frequently encountered in water bodies.[18] CBZ is efficiently taken up by many plants[5,7,8,19,20] due to its chemical properties, stability, and relative mobility in the soil.[21−23] Despite its stability in the soil, CBZ is translocated and metabolized in plants’ tissues (particularly in leaves) by plant cytochrome P450 enzymes.[5,24,25] Root uptake and reactive transport are plant- and soil-dependent,[16,23,26] making it difficult to generalize results. Numerical models play an important role in better understanding of physicochemical processes involved in the dynamic uptake of chemicals and risk assessment. To predict the behavior of non-ionic organic compounds in the soil–plant continuum, Brunetti et al.[27] coupled the widely used Richards-based solver, HYDRUS-1D,[28] with a modified version of the multicompartment dynamic plant uptake model developed by Trapp.[29] The fully coupled soil–plant model is able to provide a comprehensive mechanistic description of transport and reaction processes in the soil–plant continuum. The model was tested for three leaf vegetables (spinach, lamb’s lettuce, and arugula). The final concentrations of CBZ and its two metabolites in only two plant tissues (i.e., roots and leaves) were measured and simulated. Despite promising results, available observations were limited to draw conclusions about the model’s applicability and generalizability.

To further test the model performance, a new, more comprehensive experiment involving the exposure of green pea plants (Pisum sativum L.) to CBZ was designed and carried out in a controlled environment. The green pea plant is a widespread edible vegetable, which was chosen due to its efficiency in metabolizing CBZ[23] and its relatively fast growth and fruit formation. Extensive measurements in the soil and the plants were used in this study to calibrate the soil–plant model and assess, for the first time, its predictive uncertainty in a probabilistic Bayesian calibration framework. The calibrated model was then coupled with Global Sensitivity Analysis to identify the factors driving the accumulation of CBZ in the edible parts of the plant and shed light on the role of different physicochemical properties on the translocation of CBZ in the soil–plant continuum.

Materials and Methods

Case Study Description

Green pea plants (P. sativum L.) were planted in 24 soil columns (a height of 20 cm, a diameter of 15.4 cm) under greenhouse conditions. The soil was taken from the surface horizon (0–25 cm) of a Haplic Chernozem developed on loess (Table S1). Soil samples were first air-dried under greenhouse conditions (up to a soil–water content of 0.15 g g–1), then carefully disintegrated into aggregates with a diameter smaller than 5 mm, homogenized, and packed into the plastic columns to obtain the same bulk density (1.1 g cm–3) in all columns. Columns were then wetted to a soil–water content of approximately 0.24 cm3cm–3. Three plant seeds were sown directly into the columns in a triangular spacing. Plants (three per each column) were initially irrigated with fresh water for 16 days. Next, plants in 16 soil columns were irrigated with a solution of CBZ and in eight columns with fresh water. The columns were weighed before and after irrigation to record the water balance and estimate evapotranspiration. An evaporation pan (diameter = 19.7 cm) was used to estimate the evapotranspiration demand in the greenhouse. Irrigation doses (Table S2) depended on decreases in average soil column water contents to keep plants’ optimal conditions. No outflow at the bottom was observed. The experiment was performed between May 15 and July 2, 2018, under greenhouse conditions (natural light, air humidity of 30%–40%, and air temperature of 20–24 °C).

The columns were analyzed as follows. One column was dissembled on the 16th day to examine fresh and dry masses of stems, leaves, and roots. Next samplings were carried out on the 23rd, 30th, 41st, and 48th day. Four columns were always analyzed to obtain information about a plant growth (i.e., the mass of roots, stems, leaves, and fruits, Figure S1), concentrations of CBZ and its metabolites in plant tissues, and concentrations of compounds in four soil layers. One column irrigated with fresh water was similarly analyzed as a control. Plant leaves were scanned, and their area (Figure S3) was evaluated using the ImageJ software, version 1.47.[30] To monitor changes in the soil hydraulic properties over time, three undisturbed 100 cm3 soil samples were taken from one soil column and analyzed using the multistep outflow experiment.[31] Soil water contents, θ (at depths of 2.5 and 15 cm), and pressure heads, h (at depths of 5 and 15 cm), were measured in four columns during the entire experiment using ECH20 EC-5 sensors[32] and microtensiometers Tensior 5,[33] respectively. The Freundlich sorption isotherm and the dissipation half-live of CBZ in soil were evaluated using the batch sorption[34] and degradation[35] experiments, respectively. The resulting parameters are presented in Table .

Table 1

Model Parameters, Their Bounds, and Calculated 5%, 50%, and 95% Quantiles of the Parameters’ Posterior Distributions from the Bayesian Analysisa

			posterior distributions’ quantiles
parameterb	parameter description	parameter bounds	5%	50%	95%	S1 [-]
Soil
θr1 [cm3cm–3]	residual water content	0.18
θs1 [cm3cm–3]	saturated water content	0.51
α1 [cm–1]	VGM shape parameter	(0.001, 0.1)	0.009	0.013	0.016	0.11
n1 [-]	VGM shape parameter	(1.1, 3)	1.94	2.07	2.51	0.03
Ks1 [cm day–1]	saturated hydraulic conductivity	70
θr2 [cm3cm–3]	residual water content	0.18
θs2 [cm3cm–3]	saturated water content	0.51
α2 [cm–1]	VGM shape parameter	(0.001, 0.1)	0.015	0.022	0.035	0
n2 [−]	VGM shape parameter	(1.1, 3)	2.02	2.45	2.78	0
Ks2 [cm day–1]	saturated hydraulic conductivity	70
ρb [g cm–3]	bulk density	1.1
λL [cm]	dispersivity	1
KfCBZ [cm3β μg1-β g–1]	soil–water partition coefficient	(1.1, 10)	2.24	2.44	2.63	0.07
βCBZ [-]	Freundlich exponent	0.89
μLCBZ [day–1]	CBZ degradation rate in the liquid phase	0.0068
μSCBZ [day–1]	CBZ degradation rate in the solid phase	0.0068
P0 [cm]	Feddes’ parameter	–15
POpt [cm]	Feddes’ parameter	–30
P2H [cm]	Feddes’ parameter	–300
P2L [cm]	Feddes’ parameter	–500
P3 [cm]	Feddes’ parameter	–8000
CmaxCBZ [g cm–3]	maximum CBZ concentration taken up by roots	(1 × 10–9, 1 × 10–8)	3.75 × 10–9	4.05 × 10–9	4.38 × 10–9	0.14
CmaxEPX [g cm–3]	maximum EPX concentration taken up by roots	(0, 1.0 × 10–8)	9.0 × 10–12	8.0 × 10–11	1.6 × 10–10	0
Roots
W [cm3gfw–1]	root water content	0.88
KRW [cm3gfw–1]	roots–water partition coefficient	(1, 30)	11.8	13.3	15.1	0.03
Mmax [gfw]	maximum roots mass	308
M0 [gfw]	minimum roots mass	15
Kgr [day–1]	root growth rate	0.2
τR [day–1]	CBZ degradation rate in roots	(0, 0.55)	0.0	0.01	0.02	0.04
Stem
W [cm3gfw–1]	stem water content	0.84
KSW [cm3gfw–1]	stem–water partition coefficient	(1, 30)	10.5	11.8	12.8	0.06
Mmax [gfw]	maximum stem mass	591
M0 [gfw]	minimum stem mass	10
Kgr [day–1]	stem growth rate	0.14
τS [day–1]	CBZ degradation rate in stem	(0, 0.55)	0.04	0.05	0.07	0.06
Leaves
W [cm3gfw–1]	leaves water content	0.84
KLW [cm3gfw–1]	leaves–water partition coefficient	(1, 30)	3.12	15.2	27.3	0
Mmax [gfw]	maximum leaves mass	757
M0 [gfw]	minimum leaves mass	14
Kgr [day–1]	leaves growth rate	0.13
τL [day–1]	CBZ degradation rate in leaves	(0, 0.55)	0.38	0.44	0.50	0
SA [cm2gfw–1]	leaves specific area	70
Fruits
W [cm3gfw–1]	fruit water content	0.83
KFW [cm3gfw–1]	fruit–water partition coefficient	(1, 30)	2.7	15.2	26.5	0
Mmax [gfw]	Maximum fruit mass	720
M0 [gfw]	Minimum fruit mass	0
Kgr [day–1]	fruit growth rate	0.65
τF [day–1]	CBZ degradation rate in fruit	(0, 0.55)	0.0	0.02	0.04	0.08
SA [cm2gfw–1]	fruit specific area	50
Compounds
m(CBZ) [g mol–1]	molar mass of CBZ	236.27
m(EPX) [g mol–1]	molar mass of EPX	252.28

The last column reports the first-order sensitivity indices (S1) calculated using the RBD-FAST method to assess the influence of different factors on the accumulation of CBZ in the edible fruits.

The subscripts 1 and 2 indicate the first and second soil horizons, respectively.

Analytical Methods

Concentrations of compounds in soils and plant tissues were analyzed using the methods proposed and validated for CBZ and its four metabolites in our previous studies.[16,17,21,23,36] Methods are also briefly described in the Supporting Information. Soils (four soil layers per each column) and plant tissues were freeze-dried and grounded. Compounds remaining in 2 g of soil samples were extracted with 4 mL of extraction mixture 1 (acetonitrile/water, 1/1, v/v, acidified with 0.1% of formic acid) followed by 4 mL of mixture 2 (acetonitrile/2-propanol/H2O, 3/3/4, v/v/v, acidified with 0.1% of formic acid) in an ultrasonic bath (DT 255, Bandelin Electronic, Sonorex Digitec, Berlin, Germany) for 15 min. Compounds in 0.05 g of plant tissue samples were extracted with 1 mL of extraction mixture 1 (acetonitrile/water, 1/1, 0.1% of formic acid) by shaking at 1800 min–1 for 5 min (TissueLyser II, Quiagen, Germany). A triple-stage quadrupole mass spectrometer, Quantiva (Thermo Fisher Scientific, San Jose, CA, USA), coupled with an Accela 1250 LC pump (Thermo Fisher Scientific) and HTS XT-CTC autosampler (CTC Analytics AG, Zwingen, Switzerland), was used for the analysis of irrigation water. A hybrid quadrupole-orbital trap mass spectrometer, Q Exactive HF Hybrid Quadrupole-Orbitrap Mass Spectrometer (Thermo Fisher Scientific, USA), operated in high-resolution product scan mode (HRPS), was used instead of a triple quadrupole for more complex soil and plant extracts. A Hypersil Gold aQ column (50 mm × 2.1 mm i.d., 5 μm particle size, from Thermo Fisher Scientific San Jose, CA, USA) was used for the chromatographic separation of these target compounds. A detailed description of the instrument settings can be found in ref (37). Recovery and LOQs of compounds in different matrices are presented in Tables S4 and S5, respectively.

Modeling Theory

Model Description and Setup

The numerical model developed by Brunetti et al.[27] is used to simulate the translocation and transformation of CBZ in the soil–plant continuum. The model combines the widely used hydrological model, HYDRUS-1D,[28] with a multicompartment model of dynamic plant uptake.[38] For a thorough description of the model, please refer to ref (27).

The Richards equation describes the variably saturated water flow in the soil:where t is time [T], z is the vertical coordinate [L], θ is the volumetric water content [L3L–3], K is the unsaturated hydraulic conductivity [LT–1], h is the pressure head [L], and S is a sink term representing root water uptake [T–1]. The transport of the ith solute in the soil is described using the advection-dispersion-reaction equation, assuming that solutes can exist only in the solid and liquid phases:where C is the concentration of the ith solute in the liquid phase [ML–3], s is the concentration of the ith solute in the solid phase [MM–1], ρ is the soil density [ML–3], D is the dispersion tensor for the ith solute in water [L2T–1], q is the water flux [LT–1], r is the root solute uptake term [ML–3 T–1], and ϕ represents the reaction sink/source term [ML–3 T–1].

The soil column profile (a depth of 20 cm) is discretized in 100 finite elements refined at the top to accommodate gradients induced by the atmospheric conditions. The soil domain is divided into two horizons (0 ≤ z < −10 cm; −10 ≤ z ≤ −20 cm) to account for the effect of different root densities (Figure S2) on soil hydraulic properties, which are described by the unimodal van Genuchten–Mualem (VGM) function[39] (parameters: θ [L3L–3], θ [L3L–3], α [L–1], n [-], and K [LT–1]). The batch experiment data, presented by Schmidtová et al.,[34] indicates that the Freundlich adsorption isotherm can describe the sorption of CBZ to the solid phase (parameters: K [L3β M1-β M–1] and β [-]). Measured data show that carbamazepine 10,11-epoxide (EPX) is the only CBZ metabolite detected in the soil during the experiment (Table S7). The CBZ degradation process in the soil is thus simulated using a sequential decay chain with first-order rate coefficients, μ and μ [T–1], accounting for the transformation of CBZ in EPX in the liquid and solid phases, respectively.

The measured water content at the beginning of the experiment is used as the initial condition. No solute was present in the soil before starting the experiment. An atmospheric boundary condition is applied at the soil surface, while a seepage face boundary condition (BC) is set at the bottom. The concentration flux across the top boundary is simulated using the classic Cauchy-type BC, while a zero concentration gradient is imposed at the bottom.

The piecewise linear Feddes function parameterizes the root water stress[40] and regulates actual transpiration (parameters: P0, POpt, P2H, P2L, and P3). Simultaneously, a root density function modulates actual root water uptake depending on the root growth. A logistic growth function was fitted to the measured wet masses of roots at different depths (Figure S2) to determine the average root growth rate, which was then used to estimate the root density during the numerical simulation. The measured evapotranspiration demand is partitioned into potential evaporation and transpiration using a time series of leaf area indices, which were calculated by fitting a logistic function to the measured one-sided leaves area (Figure S3). The difference between total evapotranspiration, measured by weighing the columns, and transpiration is used to estimate evaporation from the soil surface (Figure S4). All compounds dissolved in water are passively taken up by plant roots when their concentrations are lower than maximum allowed solution concentrations, C [ML–3].[41]

Plants are conceptualized into four compartments: roots, stem, leaves, and fruits. Each compartment is assumed to have a logistic growth[38] (parameters: W [L3M–1], Mmax [M], M0 [M], K [T–1], and S [L2M–1]). The specific areas of fruits and leaves are used to partition transpiration fluxes. The CBZ solution was carefully applied on the soil surface to minimize its dispersion in the surrounding environment. Therefore, the effect of gaseous uptake and particle deposition could be neglected in the model. Similarly, volatilization is excluded from modeling due to the nonvolatile behavior of CBZ. Tissue-water partitioning (K, K, K, and K [L3M–1]) is assumed. Since the only detectable CBZ’s metabolite in plants was EPX (Table S6), the analysis was restricted to these two compounds using a first-order degradation coefficient, τ [T–1], in each compartment, which simulates the transformation of CBZ in EPX. While input concentrations and sorption parameters are converted to molar units before executing the model, output concentrations are reported in ng/g for visualization purposes.

Bayesian Analysis

This study adopts a probabilistic approach based on the Bayesian inference to calibrate the dynamic plant uptake model and assess its predictive uncertainty. Model calibration is required to estimate parameters that are difficult to measure in the laboratory (e.g., degradation rates)[3] and others (e.g., VGM parameters) that need to be adjusted to match the experimental conditions, which can be different from those encountered in the laboratory.[42] To this aim, the multimodal nested sampling (MULTINEST)[43] is combined with the soil–plant model and measured data to estimate the parameters’ posterior distribution. Multiple studies have demonstrated the accuracy and computational efficiency of this algorithm,[44−47] whose thorough description can be found in refs (43) and (48). As in refs (44) and (45), a convergence analysis is used to assess the stability and accuracy of the MULTINEST estimator.

If we assume that error residuals are uncorrelated and normally distributed with a constant variance, σ2, the log-likelihood function for the jth set of measurements can be written aswhere H(Ω) and y̅ are the ith model realization and its corresponding measured value, respectively, and k is the number of observations in the jth set of measurements. The measured data used in the Bayesian analysis include the time series of volumetric water contents (θ) and pressured heads (h) at two different locations, CBZ (C) and EPX (C) concentrations in the soil and CBZ and EPX concentrations in the roots (C, C), stem (C, C), leaves (C, C), and fruits (C, C). Therefore, the log-likelihood function L(Ω) is the aggregated sum of single likelihoods. The mean and variance of each measurement set used in the model calibration are calculated from the column replicates (Tables S6 and S7).

The dimensionality of the inverse problem is reduced by fixing specific model parameters (Table ). In particular, the results of multistep outflow experiments are used to fix the residual and saturated water contents and saturated hydraulic conductivity. This choice is motivated by the dynamics of the measured soil moisture, which never approached nearly saturated conditions during the experiment, thus providing negligible information content in that range. Feddes’ parameters for green peas (i.e., P0, POpt, P2H, P2L, and P3 in Table ) are taken from the literature.[49] Results of batch and degradation experiments are used to fix the Freundlich exponent and the first-order degradation coefficients for CBZ, respectively. The same sorption parameters are assumed for EPX, whose further degradation in the soil is neglected.[21] Water contents of plants’ tissues are assumed constant and calculated from the average measured fresh and dry weights during the experiment. Growth parameters of different plant compartments are estimated by fitting logistic functions to the time series of measured plant masses (Figures S1 and S2). Due to these prior assumptions, the dimensionality of the problem is reduced to 15 parameters. Uniform prior distributions are used in the Bayesian analysis. Parameter bounds are provided in Table .

Bioaccumulation of Carbamazepine in Edible Fruits: Global Sensitivity Analysis

After the model calibration and the uncertainty assessment, a Global Sensitivity Analysis is performed on the calibrated parameters to identify the most important factors driving the bioaccumulation of CBZ in green peas’ fruits. The analysis provides a statistical basis to shed light on the relative importance of different physicochemical processes on the translocation of CBZ toward the edible parts of the plant, which are directly linked to human and animals’ intake.

The Random Balance Designs Fourier Amplitude Sensitivity Test (RBD-FAST)[50,51] is applied. The method combines the accuracy of the classic Fourier Amplitude Sensitivity Test[52] with the computational cheapness of Satterthwaite’s random balance designs[53] to estimate the first-order effect (S1) in a variance-based context.[54] Higher S1 values are attributed to more influential parameters. The main advantage of the RBD-FAST method is that the total number of model runs is kept down to N instead of d × N (like in Sobol’ or FAST methods), where d is the number of factors investigated. In the present study, a convergence analysis is used to determine N (Figure S5).

Results and Discussion

Bayesian Analysis

Figure shows a comparison between the measured (black circles) CBZ concentrations in different plant’s compartments and soil (0 < z < −5 cm), pressure heads (z = −5 cm), and volumetric water contents at two different depths (z = −2.5 and 15 cm) and the model predictive checks (blue lines) obtained by random sampling of 100 solutions from the posterior parameter distributions. Posterior predictive checks are used to assess the effect of the uncertainty on the model predictions and to “look for systematic discrepancies between real and simulated data”,[55] thus representing a fundamental tool for model criticism. Results demonstrate that the model can reproduce water flow and solute accumulation in the soil as well as the dynamics of the translocation and metabolization of CBZ and EPX in all plant compartments, with limited uncertainty and satisfactory accuracy. The highest discrepancy is observed in the fruit and leaves compartments, for which the model overestimates the solute concentrations, though the error magnitude is low. This deviation can be attributed to dynamic variations in the plant water content and the sorption coefficient during the experiment, which are not considered in the model. The volumetric water content is well described, while the model seems to underestimate the variance in the measured pressure head. However, this tendency might be partially related to a very small measurement footprint of the tensiometers, which makes them extremely sensitive to heterogeneities in the soil’s hydraulic regime not represented in the model.

Figure 1

Comparison between the measured (black circles) CBZ (left column) and EPX (right column) concentrations in different plant’s compartments, CBZ concentration in soil (0 < z < −5 cm), pressure heads (z = −5 cm), and volumetric water contents (z = −2.5 and −15 cm), and corresponding modeled values (blue lines) obtained by random sampling of 100 solutions from the posterior parameter distributions. The error bars indicate the standard deviations of the measurements. The red line indicates the model predictions obtained by using the median solution reported in Table . The mean standard deviation used in the Bayesian analysis (σ) and the median root mean square error (RMSE) are reported in each subplot.

The calculated 5%, 50%, and 95% quantiles of the parameters’ posterior distributions are reported in Table and indicate a well-posed inverse problem characterized by low parameters’ uncertainty. The difference between the VGM shape parameters for the two soil horizons is appreciable and confirms the roots’ potential influence on water flow. In particular, smaller values of α and n calculated for the first soil layer are compatible with an increased pore heterogeneity caused by a higher root density in that part of the soil column.[56] The estimated soil–water partition coefficient, K, is in line with the results of the batch experiment (i.e., K ≈ 2.8 cm3β μg1-β g–1) and with other studies,[27,34] thus confirming the overall reliability of the analysis. The narrow confidence intervals estimated for the soil hydraulic parameters suggest that the soil behavior influences the uptake and translocation processes. This aspect is further investigated in the following sensitivity analysis.

Interestingly, the Bayesian calibration framework reveals a limited solute uptake capacity of plants (Table ). This is also confirmed by the measured data, which show a much higher persistence of CBZ in the soil than in the plants (Figure S7). In particular, the CBZ and EPX median maximum concentrations taken up by plants are 4.45 and 0.1 ng cm–3, respectively. The reduced uptake mechanism is potentially explained by the difference between root permeability to chemicals and water, making the uptake of CBZ and EPX slower than that of water. The use of a maximum uptake concentration used in HYDRUS-1D shares some conceptual and mathematical similarities with the retardation factor proposed by Gredelj et al.[57] to simulate the translocation of perfluoroalkyl acids in red chicory. However, this model parameter lacks physical meaning and theoretical background to generalize it for other modeling circumstances. The difference between C and C suggests that CBZ is taken up to a greater extent than EPX, whose occurrence in plant tissues results from in-plant metabolism rather than plant uptake. This is also confirmed by the negligible amount of EPX observed in the soil (Table S7).

Comparing the estimated first-order degradation rates further reveals that most of the CBZ metabolization occurs in the leaves. This is further highlighted by the significantly higher EPX concentrations measured in the leaves compared to other compartments. In contrast, roots and stem mainly act as reservoirs for CBZ, which is either sorbed or translocated upward. This behavior was observed in other studies as well.[2,16,17,20,23,25,27] Both the experiment and the model indicate a low accumulation of a chemically active compound, such as EPX, in the plant’s edible part. The lower accumulation of both compounds in fruits can be partially explained by the significantly shorter exposure to the contamination and lower transpiration of fruits than leaves.[23] In addition, the metabolization of CBZ in fruits has been proven to be generally less efficient.[8,19,23]

The comparison between the estimated partition coefficients confirms that sorption occurs mainly in the roots and stem, while the wide confidence intervals for K and K imply the lack of their statistical identifiability. The estimated median value for roots (i.e., K = 13.3 cm3gfw–1) agrees with Chuang et al.[1] and is mainly related to the affinity of CBZ to the plant composition[58] and its hydrophobicity.[59] In particular, different plant components such as lipids, carbohydrates, waxes, lignin, and suberin can act as sorption sites for chemicals.[60] Their distribution is tissue-dependent,[61] thus leading to different sorption coefficients in different plant compartments.

Global Sensitivity Analysis

The convergence analysis indicated that 3000 model executions were sufficient to obtain stable estimations of the first-order sensitivity indices (reported in the last column of Table ) for all investigated parameters (Figure S5). Nine factors out of 15 exhibit an appreciable influence on the translocation of CBZ toward fruits. The most influential parameter is the maximum concentration taken up by plants, which regulates the solute amount entering the roots. Once the compound enters the plant, the sorption and degradation processes in the roots and stem influence its translocation into the fruits, as suggested by the calculated sensitivity indices. Higher sorption in these compartments retards the translocation of CBZ in the transpiration streams and provides more time for the in-place enzymatic degradation. On the other hand, the leaves’ role is insignificant as the phloem flux is neglected in the numerical simulations.

The role of soil hydraulic properties holds more interest. The sensitivity analysis reveals that the VGM shape parameters of the first soil horizon and the soil–water partition coefficient strongly influence the accumulation of CBZ in the fruit compartment. In contrast, the contribution of the second soil layer is insignificant. The roots and solute distribution in the soil are both contributing factors to this behavior. Root water and solute uptake occurs mainly in the upper part of the soil profile, where the root density is higher, and CBZ is more bioavailable due to the nonlinear sorption process (Figure S6). The solute movement toward deeper soil horizons is restricted by low water flow velocities and solute dispersivity. Under such conditions, the air-entry pressure parameter α1 partially regulates actual root CBZ uptake (Figure ), the CBZ concentration in the fruit compartment, and the root zone pressure head and concentration, which increase with increasing values of α1 and K, respectively. Root solute uptake is negatively correlated with α1 [(C) in Figure ], mainly due to the change in the transpiration pattern induced by different pressure head distributions in the root zone [(E) in Figure ]. In the simulated scenarios, higher values of α1 lead to a lower (in absolute value) simulated root zone pressure head, which drops below the optimal pressure head value for transpiration (i.e., POpt in the Feddes’ model), thus reducing actual root water and solute uptake. This physical behavior is exacerbated after t = 30 days [(C) and (E) in Figure ] and results in lower translocation of CBZ toward the fruit compartment [(A) in Figure ]. Conversely, for α1 = 0.001 (1/cm), the root zone pressure head is shifted toward much higher values, where the actual transpiration rate is also reduced, though the effect is less pronounced. For α1 = 0.013 (1/cm) (i.e., median value in Table ), the root zone pressure head oscillates in the optimal transpiration range (i.e., P2H < h < POpt). Therefore, it can be concluded that soil hydraulic properties affect the accumulation of CBZ in fruits by simultaneously driving water flow in the unsaturated domain and modulating the transpiration pattern. These effects then propagate nonlinearly into the plant. However, it must be emphasized that these conclusions are not general but restricted to the modeling scenario investigated in the present study.

Figure 2

Sensitivity analysis results. (A) and (B) Simulated CBZ concentrations in the fruits as a function of α and K, respectively. (C) and (D) Simulated root solute uptake as a function of α and K, respectively. (E) Simulated average pressure head (pF = log(|h|) in the root zone as a function of α. (F) Simulated average CBZ concentrations in the liquid phase in the root zone as a function of K. The red dashed lines indicate the Feddes’ parameters (Table ). The plots were obtained by performing numerical simulations with three different values of α1 and K. In particular, the median solution reported in Table was used as a reference, and only the values of α1 and K were alternatively changed to match the minimum and the maximum listed in Table .

On the other hand, CBZ sorption in the soil affects linearly simulated CBZ concentrations in the fruit compartment. Generally, higher sorption (which for CBZ increases with increasing organic matter content in soils[34]) reduces the solute availability in the root zone [(F) in Figure ], thus reducing the amount of solute entering the roots [(F) in Figure ] and the concentration in the fruits [(B) in Figure ]. Klement et al.[23] found positive correlations between K and normalized concentration loads of the parent compound CBZ and its metabolites. In particular, the correlation was significant (95% level) for roots and stems but not significant for leaves and fruits. This suggests that organic-rich soils might reduce the translocation of CBZ and similar compounds to the edible parts of the plant.

Strengths and Limitations

We consider this study to be innovative in two main aspects:

To our knowledge, this is the first study involving numerical simulations of dynamic plant uptake of chemicals that includes extensive spatial and temporal measurements of soil- and plant-related quantities. Indeed, compared to previous studies,[62] the dataset includes measured volumetric water contents and soil pressure heads at different depths, providing further information about the unsaturated soil domain. This significant and variegated amount of data allows us to characterize water flow and reactive solute transport in the soil–plant continuum with low predictive uncertainty and satisfactory accuracy.

Compared to other inverse modeling studies,[44,63−66] this study applies, for the first time, a Bayesian calibration framework to estimate both soil and plant parameters, thus providing a statistical basis for the uncertainty assessment in the soil–plant continuum. Furthermore, using the Nested Sampling algorithm allows us to estimate the posterior parameter distribution and the marginal likelihood simultaneously. While not being of immediate interest to the present study, the latter can be used to assess model complexity, a valuable tool when dealing with highly complex and parameterized models. This methodological feature represents an advantage against other approaches focused on the use of Monte Carlo Markov Chain algorithms.[63−65]

However, these strengths are paired with several limitations. First, for most practical risk assessment applications, it is not possible to obtain comprehensive measurements due to time constraints and budget limitations. This raises some questions regarding the model’s applicability in a predictive setting in case of limited observations. There is no clear and definite answer to this. Still, we advocate the ubiquitous use of uncertainty analysis (e.g., Bayesian analysis) to identify the most justifiable model complexity level in light of the information content of the data available. Furthermore, a Global Sensitivity Analysis can, as demonstrated in this study, be used in a prognostic way to target the most important processes and focus the experimental and modeling effort on them. In such circumstances, a cross-disciplinary approach is of crucial importance.

Other model limitations include the macroscopic description of root water uptake and the lack of the physical meaning of the maximum uptake concentration. The former simplifies the root system hydraulics, which plays an important role in root water and solute uptake.[67] The latter hinders the generalization of the particular modeling approach, which is appealing from a computational point of view due to its simplicity, to other conditions. Therefore, it is important to improve the representation of these processes in the model to increase its generalizability and accuracy.

Implications and Future Outlook

This study’s main aim was to evaluate further the performance of the coupled soil–plant model developed by Brunetti et al.[27] against measurements from a comprehensive experiment on the translocation and transformation of CBZ in green pea plants in a controlled environment. A probabilistic Bayesian framework was used to calibrate the model and assess its predictive uncertainty. Results confirmed that the model could reproduce accurately and with low uncertainty the transport of CBZ in the soil and its translocation and transformation in different plant tissues. The following Global Sensitivity Analysis has further highlighted the role of soil in bioaccumulation of CBZ in the fruit compartment, mainly by regulating actual transpiration streams and CBZ bioavailability in the root zone.

Overall, the study confirms that the model can be successfully used for simulating the translocation of neutral compounds in the soil–plant continuum in a partially mechanistic way. This opens new possibilities for a more comprehensive assessment of the contamination risk, allowing the development of reliable mitigation strategies for environmental pollution problems in both the soil and plant domains. On this basis, future studies should explore the possibility of extending the model applicability to ionic compounds frequently encountered in contaminated sites. However, this would require extensive modifications to the plant model.[68]