A novel multiscale biophysical model to predict the fate of ionizable compounds in the soil-plant continuum

Soil pollution from emerging contaminants poses a significant threat to water resources management and food production. The development of numerical models to describe the reactive transport of chemicals in both soil and plant is of paramount importance to elaborate mitigation strategies. To this aim, in the present study, a multiscale biophysical model is developed to predict the fate of ionizable compound in the soil-plant continuum. The modeling framework connects a multi-organelles model to describe processes at the cell level with a semimechanistic soil-plant model, which includes the widely used Richards-based solver, HYDRUS. A Bayesian probabilistic framework is used to calibrate and assess the capability of the model in reproducing the observations from an experiment on the translocation of five pharmaceuticals in green pea plants. Results show satisfactory fitting performance and limited predictive uncertainty. The subsequent validation with the cell model indicates that the estimated soil-plant parameters preserve a physically realistic meaning, and their calibrated values are comparable with the existing literature values, thus confirming the overall reliability of the analysis. Model results further suggest that pH conditions in both soil and xylem play a crucial role in the uptake and translocation of ionizable compounds.


Introduction
The latest data indicate that soil pollution represents an alarming environmental issue (Das and Horton, 2018;Eugenio et al., 2018;The Lancet Planetary Health, 2017). The European Environmental Agency identified approximately 3 million potentially polluted sites in Europe (Van Liedekerke et al., 2014), while more than 19% of Chinese soils are categorized as contaminated (CCICED, 2015). Worrying numbers are also reported in the United States and Australia (DECA, 2010), and the lack of data from developing countries suggests an underestimation of the problem.
Sustainable practices, such as the use of reclaimed wastewater and biosolids for irrigation and fertilization, respectively, constitute preferential pathways for the introduction of PPCPs in the environment, as conventional wastewater treatment processes are only moderately effective at removing these chemicals from wastewater (Berendonk et al., 2015;Christou et al., 2017aChristou et al., , 2017bReichel et al., 2014). Antibiotics (e.g., sulfonamides and macrolides) (Aznar et al., 2014;Conde-Cid et al., 2018), analgesics (e.g., non-steroidal anti-inflammatory drugs) (Mejías et al., 2021;Qu et al., 2021), and antidepressants (e.g., serotonin reuptake inhibitors) (Malvar et al., 2021;Wu et al., 2010) are frequently detected in soils, where they can be degraded into active metabolites Kodešová et al., 2016), migrate towards groundwater (Zhang et al., 2021), or be taken up by plants (Chuang et al., 2019;González García et al., 2019;Klement et al., 2020;Sabourin et al., 2012). A better understanding of the physicochemical processes driving the fate of pharmaceuticals in the soil-plant continuum is of paramount importance to assess environmental risks and develop reliable mitigation strategies.
Numerical models play a crucial role in conceptualizing and transferring to a broad audience of decision-makers the key findings of experimental campaigns. A well-calibrated soil-plant model can be used to perform long-term scenarios to quantify multiple environmental risk factors (e.g., groundwater contamination, human risk exposure), and develop and test mitigation strategies (e.g., phytoremediation). Recently, Brunetti et al. (2019) coupled the widely-used Richards-based solver, HYDRUS (Šimůnek et al., 2016), with a modified version of the multi-compartment dynamic plant uptake model developed by Trapp (2007) to simulate the reactive transport of neutral compounds in the soil-plant continuum. The model was successfully used to predict the fate of carbamazepine in green peas under transient conditions (Brunetti et al., 2021). However, being limited to neutral compounds, the model cannot be used to analyze the fate and transport of a broad spectrum of pharmaceuticals, which exist in ionic form in the environment (Polesel et al., 2015b;Prosser et al., 2014). These compounds pose serious environmental risks. Arp et al. (2017) presented a procedure to screen and rank neutral, ionizable, and ionic organic compounds according to being sufficiently persistent and mobile (PMOCs) to survive treatment processes and enter drinking water sources. Only 9% of neutral substances received the highest PMOC ranking, compared to 30% of ionizable compounds and 44% of ionic compounds. This further underlines the importance of developing numerical tools able to simulate the fate of electrolytes in the environment.
To fill this scientific gap, the model proposed by Brunetti et al. (2019) is modified in this study to take into account processes such as chemical dissociation, electrical interaction with plant biomembranes, and phloem transport. The modeling framework connects a multi-organelles model to describe processes at the cell level with a general soil-plant model to simulate the reactive transport of electrolytes in the soil and plant. The model is calibrated using a probabilistic approach against extensive measurements from an experiment on the uptake and translocation of five pharmaceuticals (i.e., citalopram, clindamycin, fexofenadine, irbesartan, and sulfamethoxazole) in green pea plants. Selected pharmaceuticals include cations, anions, zwitterions, and ampholytes, thus being representative of a broad range of compounds frequently encountered in environmental problems (Bai et al., 2021;Hawker et al., 2013;Osorio et al., 2016;Schaffer and Licha, 2015). Furthermore, this physicochemical diversity guarantees comprehensive testing of the proposed modeling approach. The predictive performance of the model that for the first time couples a mechanistic vadose zone model with a multi-compartment plant uptake model to simulate the fate of electrolytes in the soil-plant continuum, the physical meaning of its estimated parameters, as well as the strengths and limitations of the prosed approach are discussed.

Case study description
The experimental setup was described before by Brunetti et al. (2021). Briefly, the soil material was collected from the topsoil (0-25 cm) of a Haplic Chernozem developed on loess. The soil was air-dried (a final soil water content was around 0.15 g g − 1 ) and disintegrated (a diameter of aggregates was smaller than 5 mm). The soil was then packed into 24 plastic columns with an inner diameter of 15.4 cm, and equipped at the bottom with a thin non-woven textile with a 1 cm steel mesh. The final height of the soil columns was 20 cm, and a bulk density was 1.1 g cm − 3 . The soil columns were then wetted to a soil water content of approximately 0.24 cm 3 cm − 3 . The experiment took place in a greenhouse (natural light, air humidity of 30-40%, and air temperature of 20-24 • C) between May 15 and July 2, 2018. Three seeds of green pea plants (Pisum sativum L.) were sown directly into the columns in a triangular spacing. Initially, all columns were irrigated with fresh water. After 16 days, irrigation in 16 columns was changed from freshwater to a solution with 5 ionizable pharmaceuticals (citalopram, clindamycin, fexofenadine, irbesartan, and sulfamethoxazole, Table S4). The remaining columns were further irrigated with fresh water.
The choice of using spiked freshwater is not uncommon in plant uptake studies (e.g., Goldstein et al., 2014;Gredelj et al., 2020bGredelj et al., , 2020aHurtado et al., 2016;Wu et al., 2010), and it aims to limit uncertainties stemming from other components (e.g., DOC, nutrients), which could undermine the model validation. Measured irrigation doses and concentrations are listed in Table S5. Solutions for irrigation were prepared when needed and stored at 4 • C. The intended concentration of each compound in solution was 100 ng cm − 3 . Actual measured concentrations (Table S5) were fifty to one hundred times over environmentally relevant concentrations, e.g., in wastewater presented by Golovko et al. (2014aGolovko et al. ( , 2014b. Such concentrations were used to enhance the detection and quantification of all compounds and their metabolites in all matrices. Similar or even higher concentrations were also used in previous studies (e.g., Chung et al., 2017;Klement et al., 2020;Kodešová et al., 2019a;Montemurro et al., 2017).
Irrigation rates were adjusted to ensure optimal moisture conditions for plant growth and avoid the solute leaching from the soil columns. Pressure heads, h (at depths of 5 and 15 cm), and soil water contents, θ (at depths of 2.5 and 15 cm), were measured in four columns during the entire experiment using microtensiometers Tensior 5 (METER Group Inc. USA, 2018a) and the ECH20 EC-5 sensors (METER Group Inc. USA, 2018b), respectively. The columns were also weighed before and after irrigation to record the water balance and estimate evapotranspiration. An evaporation pan was used to estimate the evapotranspiration demand in the greenhouse. No outflow at the bottom was observed.
Selected columns were sequentially analyzed on the 16th (one column), 23rd (four columns), 30th (4), 41st (4), and 48th day (4) of the experiment. The above-ground plant parts were divided into stems, leaves, and fruits. The columns (i.e., soil and roots) were cut into 4 equal layers. All plant samples were weighted to characterize plant growth (i. e., the fresh and dry masses of stems, leaves, fruits, and roots in 4 column layers, Fig. S1 and Fig. S2). In addition, plant leaves were scanned, and their area (Fig. S3) was evaluated using the ImageJ software, version 1.47 (Schneider et al., 2012). Concentrations of all compounds and their metabolites in all plant samples and all four soil layers were measured using the method described below. One column irrigated with fresh water was always similarly analyzed as a control. The remaining columns were used to measure soil hydraulic properties. However, parameters of the soil hydraulic functions (see Section 2.4.2) were additionally optimized using water flow data (irrigation, evapotranspiration, soil water contents, and pressure heads) and numerical inversion, which was in detail described by Brunetti et al. (2021).

Analytical methods
The methods to analyze concentrations of studied compounds and their metabolites (clindamycin sulfoxide, N-desmethylcitalopram, N1acetylsulfamethoxazole, N4-acetylsulfamethoxazole) in plant tissues and soils were described and validated in our previous studies Klement et al., 2020;Koba et al., 2016;Kodešová et al., 2020Kodešová et al., , 2019aKodešová et al., , 2019b. The analysis was directed towards metabolites that are commonly measured in wastewater or sewage sludge and have been measured in soils and plants (e.g., Koba et al., 2017;Kodešová et al., 2019bKodešová et al., , 2019a. Details about the extraction procedure and validation of the analytical method for carbamazepine and its metabolites were also presented by Brunetti et al. (2021). Methods are also briefly introduced in the Supplementary Material. Plant tissues and soil samples were freeze-dried and grounded. Compounds accumulated in 0.05 g of plant tissue samples were extracted with a 1 mL extraction mixture (acetonitrile/water, 1/1, 0.1% of formic acid) by shaking at 1800 min − 1 for 5 min (TissueLyser II, Quiagen, Germany). 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 with 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. Concentrations of compounds in irrigation water were analyzed using 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). Concentrations of compounds in plant and soil extracts were measured using 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). 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. More detailed information about the instrument settings was presented by Grabicova et al. (2018). Information about conditions of analysis, including gradient elution conditions, m/z values, and retention time, is provided in Table S1. Recovery and Limits of Quantification (LOQs) of compounds in different matrices are presented in Table S2 and Table S3, respectively.

Modeling theory
The problem of predicting the transport of pharmaceuticals in the soil-plant domain is tackled at two different scales and based on the use of two models. The first one, originally developed by Trapp (2000) and then further refined (Trapp, 2009(Trapp, , 2004, provides a mechanistic description of the uptake and redistribution of dissociating electrolytes in plant cells' organelles, and is used to estimate solute uptake factors and plant-water partition coefficients. These parameters are then used in a soil-plant model (Brunetti et al., 2019), simulating the reactive transport of chemicals in the soil and their translocation in the soil-plant continuum. In the following sections, we first describe the theoretical backgrounds of the two models, the differences with their original implementations, and finally their joint use in the present research study.

Cell model
To predict the uptake of ionizable compounds in plants, Trapp (2000) developed a numerical model that simulates the redistribution of chemicals among different cell's organelles. The model, which was successfully applied in other studies (e.g., Buchholz and Trapp, 2016;González García et al., 2019;Trapp et al., 2008;Trapp and Horobin, 2005), divides the plant cells into four compartments: cytoplasm, vacuoles, phloem, and xylem. To enter cells, chemicals in soil pore water have to cross a negatively charged biomembrane called plasmalemma (− 120 mV). Inside, the cytoplasm (neutral pH) surrounds vacuoles, which occupy the largest fraction of the cell volume. Vacuoles are acidic and enveloped in a positively charged biomembrane called tonoplast (+20 mV). After crossing the multiple cells and the Casparian strip, water and nutrients are redistributed to different plant tissues by specialized cells, phloem and xylem, which constitute the vascular system of the plant. Phloem is highly alkaline (pH≈8.0), while xylem is generally acidic (pH≈5.5), though its pH can vary substantially due to multiple environmental stressors (e.g., Bahrun et al., 2002;Dodd et al., 2003;Wilkinson and Davies, 1997). Conceptually, the soil pore water is an outside compartment for the cytoplasm, which, in turn, is the outside compartment for vacuoles, xylem, and phloem.
Neutral compounds in the liquid phase, which are not affected by pH conditions and electrical interactions with biomembranes, can easily reach phloem and xylem and redistribute to different plant tissues. Conversely, electrolytes can dissociate into ionic forms, which interact with cell membranes. To account for these processes, the transport of electrolytes across multiple biomembranes and their accumulation in organelles is numerically described by the Nernst-Planck equation. The total unit net flux of neutral and dissociated molecules J [M s L − 2 T − 1 ] is defined in terms of activity a [M s L − 3 ] as: where the subscripts i and o indicate the flux from inside and outside, respectively; r is the number of dissociation processes; P n and P dj are membrane permeabilities [LT − 1 ] for the neutral and jth dissociated form, respectively; a n and a dj are the activities of the neutral and jth dissociated form [M s L − 3 ], respectively, which are calculated with the Henderson-Hasselbach equation. The subscripts used with units are as follows: s represents a solute, l stands for liquid, p for protein, and fw for fresh weight. N j is defined as follows: where z j is the valence of the jth dissociated form [-] (+1 cationic, − 1 anionic), E is the biomembrane potential [mL 2 T − 3 I − 1 ], F is the Faraday constant (96,484.56 C mol − 1 ), R is the universal gas constant (8314 J mol − 1 K − 1 ), and T is the absolute temperature [K]. The biomembrane permeability for the neutral molecules P n is calculated from the compound lipophilicity (Trapp, 2009). From diffusion velocities and membrane thickness, the following equation was derived: where K OW is the octanol-water partition coefficient [-]. The same equation is used for the permeability of dissociated molecules P d , for which log K OW is 3.5 log units lowered.
While the activity a is the driving force in the diffusive exchange across biomembranes, experiments often measure the total compound concentration, which includes neutral and dissociated molecules that can be dissolved or adsorbed. In particular, sorption to lipids and proteins is considered (González García et al., 2019), leading to the following expression for the sorption coefficient K [L 3 M s − 1 ]: ] usually assumed to be 0.00122 m 3 kg l − 1 (Trapp, 2007), L is the lipid fraction [M l M fw which is used to approximate the sorption of molecules to plant proteins. The flux J can be expressed in terms of fractions f [-], which relate the activity to the dissolved concentration (Trapp et al., 2008): where fractions f are calculated by considering that the activity a is lower than the liquid concentration C [M s L − 3 ] due to the molecular interactions in non-dilute solutions, which reduce the chemical potential. The Setchenov and Debye-Hückel equations are used to calculate the activity coefficients that relate the ion concentration in the liquid phase to its activity (Trapp, 2009). The same formalism is applied to ampholytes and zwitterions (Trapp and Horobin, 2005). Finally, the mass balance equations for different cell organelles can be written as: , and J is the total net flux of the chemical across the biomembrane [M s L − 2 T − 1 ]. By assuming equilibrium, it is possible to calculate the partition coefficients between organelles K [-] and soil pore water: where C W is the chemical concentration in soil pore water [M s L − 3 ].

Soil-plant model
The numerical model developed by Brunetti et al. (2019) for neutral compounds was further modified to simulate the translocation and transformation of electrolytes in the soil-plant continuum. The model combines the widely-used hydrological model, HYDRUS (Šimůnek et al., 2016), with a multi-compartment model of dynamic plant uptake (Trapp, 2007). HYDRUS numerically solves the Richards and the advection-dispersion-reaction equations to describe water flow and reactive solute transport in variably saturated porous media, respectively. For a thorough description of the model structure, please refer to Brunetti et al. (2019). In the present study, for conciseness, we report only the key equations necessary to interpret results.
The model was originally developed for neutral compounds with the assumption of phase equilibrium between root tips and soil, which implies that the inflow concentration in roots' is equal to the concentration in the soil pore water, thus leading to the following general mass balance equation for the roots compartment: , and τ R is the first-order metabolization coefficient The root-soil equilibrium assumption does not hold for electrolytes, which can dissociate into neutral and ionic forms depending on the pH of the surrounding environment and interact with the electrically charged root cell biomembranes. As a consequence, the inflow concentration in roots xylem can be different from soil pore water concentration, and actual solute uptake must be corrected: thus leading to: where. F is a positive dimensionless solute uptake correction factor [-] that mimics the effects of the interaction between electrolytes and the root membrane, as well as reduced potential uptake due to slow permeability for polar compounds. These processes influence root solute uptake and, consequently, solute transport in the soil. Therefore, the HYDRUS source code has been modified to include the solute uptake correction factor F [Eq. (9))]. The concentration at the outflow of the xylem is in flux equilibrium to roots. Since root cells are composed of cytosol and vacuoles with density ρ cell [M fw L − 3 ] equal to 1000 kg/m 3 , the root-water partition coefficient can be calculated from Eq. (7): Mechanistically, the solute uptake correction factor F can be defined as: where K XylS is a partition coefficient between soil solution and xylem [-] [Eq. (7)]. λ is another correction factor [-] that accounts for the slow permeability of polar compounds, which is described by the following equation: where P R,water is the root permeability towards water [LT − 1 ]. In the present study, the actual transpiration rate reached approximately a peak of 0.005 m 3 /day m 2 at day 38 (Fig. S4), with a corresponding root mass of 308 g (Table S18). Meier et al. (2013) report an average root specific area of 0.086 m 2 /g for green pea plants. Therefore, the estimated roots permeability to water is approximately 2.2e-09 m/s. Nevertheless, it must be emphasized that P R,water varies significantly as it is influenced by environmental (e.g., transpiration demand) and plant specific factors (Aston and Lawlor, 1979). P R,chem is the root permeability towards chemical: where f n and f dj are the fractions of neutral and dissociated molecules [-], respectively; P n and P dj are the permeabilities to neutral and dissociated molecules, respectively [LT − 1 ], which can be calculated with the following equation: where P wall is the permeability of the cell wall towards water, which is approximately 0.25 mm/s (Trapp, 2000). The same equation is used for the permeability of dissociated molecules P d , for which log K OW is 3.5 log units lower (Trapp and Horobin, 2005). Recent studies have shown that the liposome-water partition coefficient (K lipW ) can be a better predictor than log K OW to estimate sorption of ionizable chemicals (Bittermann et al., 2014;Tanoue et al., 2017). However, its use in plant uptake modeling studies has been limited. Therefore, the present modeling framework follows the original approach proposed by Trapp (2009). Before proceeding with the analysis, the cell model has been validated against the original implementation of Trapp (2009). Details about the model validation are reported in the Supplementary Material.

Phloem transport
Electrolytes in roots' xylem are translocated upward with the transpiration stream Q Xyl first to the stem, and then to leaves and fruits according to their surface areas (Brunetti et al., 2019). Conversely, phloem transports photosynthesized nutrients and chemicals downward from leaves to the stem, and then primarily to fruits. A smaller fraction can also reach roots (González García et al., 2019). Since phloem transport is of particular importance for weak acids due to the ion trap effect (Trapp, 2009), the soil-plant model of Brunetti et al. (2019) was further modified to include this process. The phloem flux Q Phl from leaves is proportional to the change in the dry fruit mass M F and can be written as: where Q Xyl is actual transpiration [L 3 T − 1 ], φ is the fraction of actual transpiration that constitutes the phloem transport into roots [-], and T C, Ph is the phloem flux per the fruit dry mass [L 3 M dw − 1 ], which is assumed to be 10 L/kg (Trapp, 2007). Eq. (10) can be rewritten as: ], which can be calculated with the following equation: A similar formalism can be applied to describe the mass balance equations for other compartments by taking into account that no further translocation with xylem is assumed in leaves and fruits, while the phloem flux is negative for the leaves compartment and positive for the stem compartment [Eq. (16)]. The stem flux is further divided into two components leaving the stem and entering roots and fruits. Eventually, volatilization, gaseous uptake, and particle deposition processes can be embedded in the modeling framework, as described in Brunetti et al. (2019). HYDRUS can simulate intraday transpiration fluctuations if high-resolution climatic data are provided. This is generally done by partitioning the evapotranspiration demand calculated with the Penman-Monteith equation. However, the plant model assumes that the phloem transport depends on the fruit mass and is always a fraction of actual transpiration [Eq. (16)]. This ratio remains constant, even if HYDRUS is able to simulate intraday fluctuations. Other more precise numerical approaches are available to simulate an increase in the phloem flux at night (e.g., Hölttä et al., 2006), but they increase the overall model complexity. A potential solution is to adopt differentiated xylem/phloem ratios between day and night in Eq. (16). However, high temporal resolution data are needed to test this approach.

Metabolization
The soil-plant model considers compound degradation in both soil and plants compartments. A comprehensive description is provided in Brunetti et al. (2019). Here, we provide only key equations. HYDRUS can simulate sequential degradation chains from parent to daughter species using first-order rate coefficients µ [T − 1 ] in the liquid and solid phases. The plant metabolization is modeled using one metabolization matrix Γ j for each compartment, which consist of multiple first-order rate coefficients providing connections between parent and daughter species. The matrix Γ j is a square matrix of a dimension N (i.e., a number of metabolites considered): If only the parent compound is modeled, then Γ j reduces to a single degradation coefficients, which describes the transformation of the parent compound into unknown metabolites or its mineralization.

Model calibration and uncertainty assessment
Five inverse modeling scenarios were used to assess the performance of the model in reproducing experimental data involving uptake and translocation of five pharmaceuticals [i.e., citalopram (CIT), clindamycin (CLI), fexofenadine (FEX), irbesartan (IRB), and sulfamethoxazole (SUL)]. Clindamycin sulfoxide (CLIS), a metabolite of clindamycin, was also quantified in the soil and plant samples. Therefore, this metabolite was included in the inverse modeling scenario CLI. Concentrations of metabolites of other pharmaceuticals were mostly below the limits of quantification. Therefore, they are not considered in this study, and the analysis is restricted to the parent compound. In the following sections, details about the calibration strategy, model setup, and algorithms are provided.

Calibration strategy
Theoretically, the previously described cell and soil-plant models can be solved sequentially. The former model is used first to calculate the uptake correction factor F and the partition coefficients between plant compartments and xylem (K PXyl ) and phloem (K PPhl ). This information is then passed to the latter model to estimate translocation and redistribution of electrolytes in different plant compartments. However, this approach would formally require calibrating multiple parameters in the cell model, thus significantly increasing the dimensionality of the inverse problem.
To circumvent this problem, we first calibrate the solute uptake correction factor and the partition coefficients in the soil-plant model for multiple pharmaceuticals. Then, we validate and discuss the meaning and implications of these partition coefficients using the cell model. In this way, we reduce the dimensionality of the calibration problem and the associated computational burden while still retaining the mechanistic traits of the cell model in a validation step.

Model description and setup
The soil domain discretization in HYDRUS, as well as the boundary and initial conditions, for the same experimental setup are described in Brunetti et al. (2021) and not repeated here for conciseness. As in Brunetti et al. (2021), soil hydraulic properties are described by the unimodal van Genuchten-Mualem (VGM) function (van Genuchten, 1980) The piecewise linear Feddes function parameterizes the root water stress (Feddes et al., 1978) and regulates actual transpiration (parameters: P0, POpt, P2H, P2L, P3). Simultaneously, a root density function modulates actual root water uptake depending on the root growth, which is assumed to be logistic. A time series of Leaf Area Indices (LAI) is obtained by fitting a logistic function to the measured one-sided leaves area (Fig. S3). LAIs are then used to partition evapotranspiration into potential evaporation and transpiration. The difference between total evapotranspiration, measured by weighing the columns, and transpiration is used to estimate evaporation from the soil surface (Fig. S4).
The batch experiment data of Kodešová et al. (2020) and Schmidtová et al. (2020) indicate that the Freundlich isotherm (parameters: (K f [L 3β M 1-β M − 1 ], β CBZ [-]) well interprets the sorption of pharmaceuticals to the solid phase, and thus this parameterization is used in the model. Sorption of particular compounds can decrease due to their competition for the same sorption sites or increase due to their cooperation. This behavior was observed by Schmidtová et al. (2020) for the five compounds analyzed in the present study. Due to an ambiguous impact of a compound mixture on sorption behavior, sorption parameters were calibrated. The compound degradation in the soil is described by the first-order degradation coefficient, µ S [T − 1 ]. The only detected metabolite in the soil during the experiment was CLIS. The transformation from parent to child compound and its subsequent transport and degradation is thus considered only for CLI.
Each plant compartment is assumed to have a logistic growth (parameters: The specific areas of fruits and leaves are used to partition transpiration fluxes. All pharmaceuticals considered in the present study have very low air-water partition coefficients, and thus volatilization, as well as gaseous uptake and particle deposition, are neglected. As for the soil, the metabolization chain in plants is considered only for CLI, whose main metabolite (CLIS) can be taken up from the soil, generated in plants' tissues, and transported in xylem and phloem.

Bayesian analysis
A probabilistic approach based on the Bayesian inference is used to calibrate the soil-plant model and assess the parameters' uncertainty. A nested sampling algorithm (Feroz et al., 2009) is combined with the model and observations to estimate the parameters' posterior distributions. A thorough description of the algorithm, which was successfully applied in multiple studies (e.g., Brunetti et al., 2020bBrunetti et al., , 2020aElsheikh et al., 2013;Liu et al., 2016), is provided in Skilling (2006) and Feroz et al. (2009). As in Brunetti et al. (2020b) and Liu et al. (2016), a convergence analysis is used to assess the stability and accuracy of the nested sampling estimator.
We assume that measurement errors are normally distributed with an average constant variance, σ 2 , and thus the log-likelihood function ℓ(Ω) for the jth set of measurements can be written as where H i (Ω) and ỹ i are the ith model realization and its corresponding observation, respectively, and k is the number of observations in the jth set of measurements. The measured data used in the Bayesian analysis include the compound concentration in the soil at four different depths and in four plant compartments (roots, stem, leaves, and fruits). 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 (see the Supplementary Material): The dimensionality of the inverse problem is reduced by fixing specific model parameters (Table S18). In particular, Brunetti et al. (2021) calibrated VGM parameters for the same experimental setup. Therefore, these parameters are set to their median values reported in Brunetti et al. (2021) and excluded from the model calibration. The plant growth parameters are obtained by fitting a logistic function to measured masses (Fig. S1), while water contents of plants' tissues are assumed to be constant and calculated from the average measured fresh and dry weights during the experiment (Fig. S1). The longitudinal dispersivity and Feddes' parameters are fixed according to Brunetti et al. (2021). Consequently, only solute uptake correction factors, and sorption and reaction parameters in the soil-plant domain, are inversely estimated.
The number of calibrated parameters is further reduced by making additional theoretical considerations. In particular, since volatilization is neglected, and no further translocation with xylem is assumed in leaves and fruits, their partition coefficients do not affect the accumulation of electrolytes in these compartments and can be set to any realistic value in their parameter space. Similarly, since the roots and fruits are terminal compartments for the phloem transport and roots exudation is neglected, their phloem partition coefficients do not influence the transport of pharmaceuticals in the plant. Consequently, the number of parameters to be calibrated is 12 for scenarios CIT, FEX, IRB, and SUL without metabolites. It is 18 for CLI, which includes an additional uptake correction factor and five additional reaction parameters for its metabolite CLIS. Due to only minor changes in the molecular structure of CLI and CLIS, their sorption coefficients are assumed to be the same.
Uniform prior distributions are used in the Bayesian analysis. The bounds of prior distributions (Table 1) for plant uptake and sorption parameters (i.e., F, K PXyl , K PPhl ) were set using the previously described cell model, whose parameters are detailed in Table S19. On the other hand, a literature survey was used to set the bounds of degradation coefficients and soil sorption parameters (Brunetti et al., 2021;Huynh and Reinhold, 2019;Kodešová et al., 2020;Schmidtová et al., 2020).

Bioaccumulation of pharmaceuticals
Before discussing the outcome of the numerical analysis, which is the main aim of this study, we briefly explain the differences in the observed bioaccumulation of pharmaceuticals in the soil-plant continuum. The role of different mechanisms on the uptake and translocation processes is further investigated in the next section about modeling results, where we combine experimental and numerical findings.
Experimental data used in the numerical analysis are shown in Fig. 1  (gray circles) and reported in the Supplementary Material (Tables S6-S17). Further, Fig. S6 in the Supplementary Material also shows measured average fractions of the applied solute mass in the soil and different plant compartments. Pharmaceuticals show important bioaccumulation differences. CIT has the highest occurrence in the soil and is efficiently translocated into plants, where it is detected at meaningful levels in leaves. This agrees with observations of Kodešová et al. (2019), who reported an appreciable translocation of CIT in shoots. This behavior can be explained by CIT's largest persistency in soil, and moderate lipophilicity (log K OW = 3.5), and poor affinity to proteins (Lin et al., 2014), which facilitate its translocation with the transpiration stream. A similar trend was observed in other studies for other selective serotonin reuptake inhibitors (Redshaw et al., 2008;Reichl et al., 2018).
In contrast, FEX and IRB were not detected in shoots, only in soil and roots. In comparison, IRB concentrations in the roots were more than twice those of other compounds and lower than CIT and FEX in soil, suggesting a strong affinity to roots' components and poor sorption to the soil matrix. Its high lipophilicity (log K OW = 6) and protein affinity (Yang et al., 2012) can partially explain its occurrence in roots and poor mobility in plants (Dettenmaier et al., 2009). Nevertheless, the presence of untracked metabolites generated by in-plant enzymatic biodegradation cannot be excluded a priori by comparing observations (Villette et al., 2019). The combined effect of cation-bridging (Polesel et al., 2015a), surface complexation (Klement et al., 2018), and hydrophobicity (log K OW = 5.6) also explain sorption of the zwitterion FEX to roots and its reduced translocation in shoots, which was observed in Kodešová et al. (2019a) and reported for similar antihistamine in other studies (e.g., Wu et al., 2012).
The macrolide CLI and its main metabolite CLIS, which was already detected in irrigation water (Table S5), exhibit the lowest accumulation in roots, while they both persist in soil. Their concentration in shoots is higher than FEX and IRB, but lower than CIT and SUL, thus suggesting a potential role of in-plant metabolism, a process that was observed for other antibiotics and macrolides in vegetables (e.g., Tian et al., 2019). Conversely, SUL displays the highest and lowest concentrations in leaves and soil, respectively. Its translocation towards shoots was already reported by other studies (Chitescu et al., 2013;Herklotz et al., 2010;Klement et al., 2020;Kodešová et al., 2019b;Tanoue et al., 2012), and can be mainly due to its low lipophilicity (log K OW = 0.89), while its negligible occurrence in soil can be a combination of efficient root uptake and degradation processes (Koba et al., 2017;Liu et al., 2010). Concentrations of other metabolites were mostly below the limits of quantification, which can be explained by very high transformation (formation and dissipation) rates (N1-acetylsulfamethoxazole, N4-acetylsulfamethoxazole) or slow transformation rates of CIT. In some cases, it could be due to larger values of the limits of quantification (Table S3). Fig. 1 shows a comparison between the measured (gray circles) and modeled average concentrations of pharmaceuticals in soil (orange lines) and in different plant compartments (green lines) obtained by random sampling of 100 solutions from the posterior parameter distributions. Results demonstrate that the model is able, except in few cases, to reproduce with limited uncertainty and acceptable accuracy the reactive transport of multiple pharmaceuticals in the soil-plant domain. Generally, model predictive checks don't indicate systematic discrepancies between model predictions and observations (Gelman et al., 2004), suggesting that the model should not be rejected and is well suited for the numerical analysis.

Model predictive performance
Root water and solute uptake occurs mainly in the upper part of the soil profile, where the root density is higher, and compounds are more bioavailable due to the non-linear sorption process (Fig. S5). The chemicals accumulation of all pharmaceuticals in soil and roots, as well as the transformation of CLI into CLIS and its subsequent translocation, are well predicted. Similar good predictive performance is observed in leaves, and partially in the stem and fruits. The model shows a high predictive uncertainty in the stem compartment for FEX, thus simultaneously suggesting that the information content of measurements is not sufficient to constrain model predictions and that the FEX concentrations in the stem are overly sensitive to small variations in roots and soil. A comparable but less evident behavior is observed for SUL in fruits, for which the model tends to predict slight compound bioaccumulation not detected in measurements. Simplifications in the definition of the phloem transport, and inaccuracies in non-calibrated model parameters (Table S18) and observations can partially explain these deviations. In particular, the quantification limit for SUL in fruits (0.19-0.63 ng g − 1 , Table S3) was one order of magnitude higher than for other compounds. In addition, LOQs for SUL metabolites are also relatively high (N1 -0.63-2 ng g − 1 , N4 -0.2-0.69 ng g − 1 ). It can be hypothesized that if LOQs for SUL and its metabolites were lower, their concentrations in fruits could be, in some cases, quantified. However, despite these problems and considering the complexity of predicting the fate of multiple electrolytes, the model exhibits satisfactory performance and guarantees a comprehensive description of reactive transport processes in both soil and plant tissues ( Fig. 1 and Fig. S5).

Parameters assessment
The calculated 5%, 50%, and 95% quantiles of the parameters' posterior distributions for all pharmaceuticals are reported in Table 1. Posterior distributions show a mixed behavior, with some parameters precisely inferred and others exhibiting an appreciable uncertainty. In the following sections, we compare estimated parameters with the literature values, discuss their physical meaning and consequences for multiple processes, and attempt to provide a mechanistic explanation by using the Cell model.

Table 1
Calibrated model parameters, their bounds, and 5%, 50%, and 95% quantiles of the parameters' posterior distributions calculated from the Bayesian analysis. CIT, CLI, FEX, IRB, and SUL refer to columns with five pharmaceuticals, i.e., citalopram, clindamycin, fexofenadine, irbesartan, and sulfamethoxazole. a CLI and CLIS are assumed to have the same sorption parameters in soil and plants but different metabolization rates.

Root solute uptake
The low uncertainty of the estimated solute uptake correction factors F indicates that the interaction between electrolytes and roots membrane plays a fundamental role in the chemical uptake and transport (Haynes, 1980;Trapp, 2009). As previously described [Eqs. (12), (7), (5)], the electrolytes fluxes across the root biomembrane are strongly influenced by soil and xylem pH, which govern the chemicals' speciation inside and outside roots, and thus their ability to permeate the root membrane. Recently, Delli Compagni et al. (2020) used a combination of Monte Carlo and Sensitivity analyses to identify factors influencing bioaccumulation of multiple compounds in plants. The analysis revealed that xylem and phloem pH had a large impact on the simulated concentrations in plants' tissues, thus confirming the findings of our present study. Fig. 2 shows the color maps obtained using the Cell model, illustrating the mutual influence of soil and root xylem pH on the uptake correction factor, F, for each pharmaceutical under consideration. The white bands correspond to the estimated 5% and 95% quantiles of the solute uptake correction factor reported in Table 1 and identify the plausible range of soil and root xylem pH values.
The narrow bands for CLI and SUL suggest that the estimated solute uptake factors can only be explained by similar pH conditions in soil and root xylem, whereas the existence domain is broader for other compounds. While commonly assumed acidic (Trapp, 2009), root xylem pH can vary considerably due to multiple environmental stressors. In particular, under drought conditions, xylem pH can significantly increase (Bahrun et al., 2002;Dodd et al., 2003;Gloser et al., 2016;Wilkinson and Davies, 1997), activating the guard cell abscisic acid receptors in the leaf epidermis and leading to stomatal closure to enable the plant to retain water (Wilkinson, 1999). However, as shown in Brunetti et al. (2021), the soil water content during the experiment oscillated between 0.20 and 0.30, and the plant growth was sustained. Therefore, root water stress was not present during the experiment and cannot be the cause of alkaline xylem conditions. Consequently, the estimated solute uptake factors in Fig. 2 can only be justified by slightly acidic conditions in both xylem and soil. Unfortunately, no soil pH measurements were performed during the experiment to confirm this conclusion. However, variations of soil conditions due to multiple stressors are not uncommon during experiments (Carter et al., 2016). Therefore, a shift from original alkaline soil conditions (i.e., soil pH ≈  (Tables S6-S17). The y-axis scale is intentionally kept the same, when possible, for multiple compounds to emphasize their different accumulation in plant tissues. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 8.0) to acidic ones during the experiment is plausible and suggested by model calibration results for multiple compounds.

Sorption and translocation in plant
Xylem pH affects not only root solute uptake but also electrolytes sorption and translocation in the plant [Eq. (11)]. Fig. 3 shows the influence of xylem pH on the root-xylem partition coefficient K RXyl (in a logarithmic scale) for each pharmaceutical. Curves (black lines in Fig. 3) are obtained using the cell model and by varying the xylem pH in the interval [4,8], while red bands correspond to the estimated 5% and 95% quantiles reported in Table 1, which indicate that partition coefficients in roots and stem are generally affected by limited uncertainty. The overlapping between curves and uncertainty bands inferred from the Bayesian analysis identifies a plausible range of xylem pH during the experiment.
Interestingly, the overlapping between curves and uncertainty bands in Fig. 3 for all pharmaceuticals exclude alkaline conditions in the xylem, and deductively in soil, thus confirming earlier discussion. In particular, while for FEX, IRB, and CIT, the plausible ranges of xylem pH are broad and extend from acid to slightly alkaline conditions, the domains of existence for CLI, SUL, and CLIS are restricted to acidic pH. This is mainly due to the molecules' fluxes between cytosol and xylem. In acidic xylem conditions, SUL is partially present in its neutral form, which can diffuse back in the alkaline cytosol, where it dissociates into anionic molecules. These molecules cannot permeate the membrane again and thus remain trapped in the cytosol. Prevalent cationic

Fig. 2.
Color maps showing the mutual influence of soil and root xylem pH on the uptake correction factor, F, for each pharmaceutical. The white bands correspond to the estimated 5% and 95% quantiles reported in Table 1. Fig. 3. Influence of xylem pH on the root-xylem partition coefficient K PXyl (in a logarithmic scale) for each pharmaceutical. The red bands correspond to the estimated 5% and 95% quantiles reported in Table 1. The y-scale is intentionally kept the same for multiple compounds to emphasize sorption differences. electrolytes, such as CLI and CIT, exhibit the opposite behavior. FEX (zwitterion) and IRB (ampholyte) sorption is characterized by a combination of cation-bridging (Polesel et al., 2015a), surface complexation (Klement et al., 2018), and hydrophobicity effects, and therefore, once in the cell, they are readily sorbed to lipids (and other components) in the vacuole. Theoretically zwitterions can form layer onto surface due to their interconnections and thus largely increase sorption (Carrasquillo et al., 2008;Vasudevan et al., 2009). Table 1 further indicates a sorption decrease in the stem, and this can be explained by a lower lipid and protein content in this compartment [Eq. (4)]. Finally, the phloem partition coefficients exhibit high uncertainty for all pharmaceuticals. There are different potential explanations for this behavior. On the one hand, the high uncertainty suggests that phloem transport played a minor role in the chemicals' translocation in plants. On the other hand, results can also indicate that observations are not informative of this process and that measurements should be enriched.
It must be emphasized that the present model, and other existing soilplant uptake models, do not explicitly account for effects such as competitive (or cooperative) sorption, which may occur when multiple compounds are simultaneously present. The main reason is that the model complexity significantly increases when a mechanistic description of such processes is considered, which may undermine the model predictive capabilities (e.g., Brunetti et al., 2020bBrunetti et al., , 2020a if extensive observations are unavailable. However, the estimated posterior distributions of sorption coefficients (Table 1), inferred from observations implicitly account for these processes, and calibrated parameters should be formally regarded as effective. This is frequent in calibration modeling studies, where the appropriate trade-off between model complexity and generalizability is difficult to identify.

Degradation in plant
The posterior distributions of degradation parameters in plant tissues show a mixed behavior (Table 1). Confidence intervals in roots are narrow for all pharmaceuticals, thus indicating their influence on the electrolytes translocation in plants. In particular, CIT and CLIS exhibit the highest degradation rates, followed by IRB, FEX, CLI, and SUL. Reichl et al. (2018) suggested that hydroxylation, demethylation, and conjugation with amino acids, were responsible for the in-plant metabolization of sertraline, a pharmaceutical belonging to the same family of serotonin reuptake inhibitors as CIT. Similarly, efficient metabolization of CLI was reported in Tadić et al. (2019). In another study, the macrolide clarithromycin, whose molecular characteristics are comparable to CLI and CLIS, was efficiently metabolized in lettuce (Tian et al., 2019). To our knowledge, there have been no studies about the metabolization of FEX and IRB in plants. However, the IRB metabolization pathway in human cells is driven by the cytochrome P450 enzymes (Stearns et al., 1995), which have also been identified in plants (Carter et al., 2018). This supports its biotransformation in plants. Conversely, FEX undergoes minimal metabolism in human cells, and therefore, its transformation in roots is not completely enzymatic. It can be speculated that the synergistic effect of roots endophytic bacteria enhanced the degradation process (e.g., Nguyen et al., 2019). In support of this, Kodešová et al. (2020) showed that the dissipation of FEX was highly correlated with the bacterial population in the soil. However, more data are needed to confirm these results. The model indicates that SUL had the lowest, but still significant, degradation in roots. This agrees with what was reported in other studies (e.g., Huynh and Reinhold, 2019;Tanoue et al., 2012;Tian et al., 2019). Li et al. (2018) showed an intensive SUL metabolism in radish tissue enzyme extracts (the SUL recoveries in root and leaf enzyme extracts after 96 h were 34% and 72%, respectively). Similarly, Wu et al. (2016) documented an even faster metabolism of SUL in carrot cell cultures.
Similar degradation rates were estimated for the stem compartment, except for FEX and CLIS, for which the confidence intervals are large, thus hindering any meaningful conclusion. In the same way, high uncertainties were calculated for IRB and FEX in leaves and all compounds in fruits. IRB and FEX do not translocate to shoots due to their high lipophilicity, and thus no data are available to infer model parameters. Conversely, degradation rates in leaves for CIT, CLI, SUL, and CLIS were well identified and generally low. This can be partially explained by low substrate levels, which cannot stimulate the enzymatic (e.g., Michaelis-Menten) biodegradation (Hurtado et al., 2016).

Soil sorption
The divergence between priors and posteriors for all compounds (Table 1) indicates that the measurements are informative for soil sorption parameters (i.e., K f and β), which consequently played an important role in the electrolytes transport and uptake during the experiment (e.g., Brunetti et al., 2021). Estimated K f and β values indicate that CIT, CLI, FEX, and CLIS strongly sorb to the soil, while IRB and SUL are mostly present in the liquid phase. This difference is mainly due to the dissociation of electrolytes (Table S4) in the soil during the experiment, which the previous discussion on plant uptake and sorption suggests is slightly acidic. In these conditions, CIT and CLI are present mainly in the cationic forms, which are attracted to the soil matrix due to its negative charge and the hydrophobicity of the compound (Schaffer et al., 2012b(Schaffer et al., , 2012a. This behavior was already discussed by Kodešová et al. (2020) and Schmidtová et al. (2020) for the same soil under laboratory conditions. The difference between the values estimated in these two studies and those reported here can be explained by the different operating conditions, pH shifts, and the influence of roots on the compound bioavailability. Conversely, SUL is present partially in anionic and neutral forms, which can explain its moderate affinity to the solid phase Schmidtová et al., 2020;Schübl et al., 2021). This was already observed in Archundia et al. (2019), who reported an increase in SUL sorption in acidic soil conditions due to specific hydrophobic interactions with soil organic carbon (Boxall et al., 2002;Kurwadkar et al., 2011;Srinivasan and Sarmah, 2014;Thiele-Bruhn and Aust, 2004). Finally, estimated IRB and FEX sorption parameters agree with what was reported in Klement et al. (2018), who indicated that the molecular structure of these two compounds plays a key role in their sorption capabilities. In particular, Polesel et al. (2015a) suggests cation-bridging effect as the driving factor in zwitterion and ampholytes sorption.

Degradation in soil
The analysis precisely infers the first-order degradation coefficients in the soil. These coefficients generally have low values, except for SUL and the metabolite CLIS, with the former exhibiting a strong dissipation, a process that was already observed in other studies (e.g., Andriamalala et al., 2018;Liu et al., 2010). Results generally agree with what Kodešová et al. (2020) reported for the same soil, with CIT exhibiting the lowest degradation followed by IRB, FEX, and CLI. Simultaneously, the model predicts that the metabolite CLIS is rapidly degraded, which contradicts what Koba et al. (2017) found in laboratory conditions. The presence of roots (Jin et al., 2019) and the different experimental conditions can partially explain this difference.

Conclusions
The main aim of this study was to develop and assess the predictive performance of a soil-plant model for the analysis of the fate and transport of electrolytes in the soil-plant continuum. The modeling framework connects a multi-organelles model describing processes at the cell level with a general soil-plant model simulating the reactive transport of electrolytes in the soil and plant. To circumvent the problem of high-dimensional model calibration, the soil-plant model was calibrated independently, and estimated parameters were then validated using the cell model.
Results are encouraging and show that the soil-plant model can predict the fate of multiple pharmaceuticals in the soil-plant system with satisfactory accuracy and acceptable uncertainty. The subsequent validation with the cell model indicates that the soil-plant parameters preserve a physically realistic meaning, and their estimated values are comparable with the existing literature values, thus confirming the overall reliability of the analysis. Model results further indicate that pH conditions in both soil and xylem play a crucial role in the uptake and translocation of ionizable pharmaceuticals. Therefore, accurate monitoring of at least soil pH is recommended for experiments focused on electrolytes transport in the soil-plant continuum to understand how dissociation processes can affect chemicals' sorption, uptake, and translocation.
The model can be further improved by refining the description of certain processes, such as plant hydraulics, changes in soil pH, sorption kinetics, etc. Experiments with spiked wastewater (or groundwater) can be used to assess the influence of other components (e.g., nutrients, DOC) on the uptake process and improve the modeling framework. However, this additional theoretical complexity must be accompanied by more comprehensive measurements (e.g., leaf water potential, metabolite tracking, etc.) to avoid degradation of the model's predictive performance. Conversely, ad hoc laboratory experiments are desirable to better characterize physicochemical processes and improve their description in the model for explanatory purposes. Nevertheless, the developed modeling framework represents a step forward in this direction and promises to be a valuable tool for practitioners and scientists assessing environmental risks.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.