Mathematical Modeling of Pilot Scale Olive Mill Wastewater Phytoremediation Units

: A mechanistic state–space model has been developed to describe the dynamics of olive mill wastewater (OMW) treatment in phytoremediation pilot units with P. granatum L. and M. communis L. plants and to assess further the relative contribution of the plants in the overall OMW remediation process. Both phytoremediation and bioremediation processes have been considered in the model, i.e., phytodegradation, rhizodegradation, accumulation of hardly biodegradable organic matter on the root tissue of plants, microbial growth, maintenance and decay, and enzymatic decomposition of organics. Maximum speciﬁc microbial growth rates for bacteria and fungi were estimated within the range of 0.164–0.236 1/h. The speciﬁc rate for the decomposition of hardly biodegradable organics both by bacteria and fungi was within the range of 10.75–72.73 mg-substrate/g-biomass · h, whereas, particularly for the high-molecular-weight polyphenols, it was 1.02–18.25 mg-substrate/g-biomass · h. The values of the transpiration stream concentration factor were greater than 0.95 for both the non-phenolic and phenolic organics, which indicates almost passive uptake of OMW organics’ mixture by the plants. The corresponding factors for inorganic N and P were estimated as greater than unity, indicating active uptake. Overall, the model predicts the experimental data well when the organic concentration of OMW is high, and it predicts that phytoremediation processes contribute by more than 91% to the removal of OMW organics and nutrients, irrespective of the wastewater organic strength.


Introduction
Olive mill wastewater (OMW) constitutes a complex wastewater with varying composition and high concentration of organics and nutrients.Recently, Mediterranean countries have individually permitted, through legislation, the direct disposal of non-treated OMW on soil [1] since no affordable and efficient treatment technology has been found to date for this type of wastewater.However, this practice imposes the risk of groundwater contamination and, at the same time, requires the use of vast land area.To eliminate the aforementioned drawbacks, Santori and Cicalini (2002) [2] introduced a phytoremediation scheme for OMW treatment by recirculation, where OMW was used as an irrigation source for planted soil in restricted areas.This technology was implemented successfully on a large scale in Italy, and the efficacy of tree-based phytoremediation for OMW treatment was well proven [3].The main advantage of this technology against physicochemical and biological methods for OMW treatment is that it can be easily implemented (without technical expertise) in small decentralized olive mills with only a minor fixed cost for soil excavation, purchase of pumping and piping equipment, and the plants to be cultivated.Petoussi and Kalogerakis (2022) [4] recently proposed the use of high economic importance of plants in OMW phytoremediation to improve the circular economy character of the previous technology.They quantified the efficiency of the proposed technology in terms of organic and nutrient removal rates and compared those rates with the corresponding removal rates of the relevant OMW-treating CWs.They also demonstrated, at the same time, the high contribution of the plants in the OMW treatment scheme.In the proposed treatment scheme, OMW constituents are removed both through processes related to bulk soil biomass activity, as well as through plant-related processes, such as plant uptake, rhizodegradation, and accumulation on plants' roots.However, the distinct contribution of each of the above processes to the overall remediation rate has not been quantified since there is no available model to describe this specific system.
To the best of our knowledge, no modeling work has been devoted yet to the description of the dynamics of OMW degradation on soil, although this practice is very common in most Mediterranean countries.In contrast, numerous studies have described thus far the simple kinetics of OMW remediation in small-scale bioreactors or sand columns [5][6][7][8][9].The most recent literature on modeling related to OMW remediation includes studies on polyphenols recovery, whereby adsorption kinetics on various matrix materials are studied both experimentally and kinetically [10][11][12].In addition, OMW flocculation-sedimentation and biological oxidation has been studied kinetically by Vuppala et al. (2019) [13], as well as co-metabolism of OMW in sequencing batch reactor under aerobic conditions after Fentonbased oxidation [14].More complex integrated mathematical models on OMW remediation are scarce and have been so far developed for processes such as OMW composting [15,16] or for the incorporation of OMW polyphenols degradation in activated sludge systems [17].
The aims of the present study are to develop a mathematical model to adequately describe OMW degradation in a pilot phytoremediation system which uses tree species with high-value products and to determine the kinetics and the relative contribution of the various bioremediation and phytoremediation processes of the proposed system in the overall remediation of OMW.Moreover, since there are many studies for direct OMW application on bare soil, where OMW is applied in doses and gradually biodegraded, it is of interest to examine whether the presence of fruit-bearing deciduous shrubs provides an advantage in terms of the overall OMW degradation rates in addition to the economic advantage of cultivating attractively priced fruits (pomegranates or myrtle berries).Furthermore, mathematical modeling is an essential tool for the optimization of the pilot system's operation and the estimation of the contribution of other important variables of phytoremediation, such as the plants' root system development, concentration of rhizosphere microbial biomass, diversity of the plant species, and wastewater synthesis.
Innovative points in the model are the conceptual division of OMW organic mixture in two pools (labile and recalcitrant), whereby microbial biomass grows on both organic fractions following Monod as well as Haldane kinetics.Moreover, OMW polyphenol concentration has been considered as a separate state variable due to the profound importance of phenolics to the recalcitrant character of OMW.

Operation of the Phytoremediation Pilot Units
The model of the present study aimed to describe the operation of the phytoremediation system presented by Petoussi and Kalogerakis (2022) [4].In the previous work, OMW was gradually remediated by recirculation through the soil in two sandy soil pilot units, each planted with two plants of either Punica granatum L. (pomegranate tree) or Myrtus communis L. (myrtle tree) species (Figure 1).Experimental data used for model calibration and validation were produced during four OMW treatment cycles performed during the summer months in two consecutive years.At the onset of each experimental cycle, both pilots were supplied with diluted OMW in the storage compartment, which subsequently was pumped continuously to the surface of the units for 12h/d and infiltrated through the soil.Raw OMW had initial average chemical oxygen demand (COD) concentration for Cycles #1-#3 (1 st year) and Cycle #4 (2 nd year), corresponding to 2700, 6200, 12,400, and 39,300 mg/L, respectively.The dynamics of the most important OMW physicochemical variables, i.e., COD, biochemical oxygen demand (BOD5), total polyphenols (TPh), total nitrogen (TN), and total phosphorus (TP), were investigated.

Model Development of the Phytoremediation Pilot Units
Due to the macroscopic character of the model, several assumptions were made to simplify hydrodynamics and mass transfer phenomena in the pilots.More specifically, the soil layer of the unit was conceptually divided into four equal-volume consecutive compartments (Figure 1b) in a tank-in-series approach [18] to account roughly for the plug-flow type of movement of the wastewater liquid phase through the soil.Assumptions on the characteristics of the plants' root system in the two pilots were based on the most relevant available information found in the literature and are described in detail in Appendix A.1.
OMW was provided to the surface of each unit by drip irrigation, where in sandy soil, the vertical water flow was reported to be much faster than the horizontal flow due to the high infiltration rate.Therefore, OMW was considered to be infiltrated only from a specific fraction of each unit's total soil (wetted soil region), shown schematically in Figure 1b.Numerical estimation of the mass of the total wetted soil region was described in detail by Petoussi and Kalogerakis (2022) [4], and its corresponding percent distribution in the unit was considered to be 15%, 25%, 30%, and 30% for the 1st to 4thcompartments, respectively.Wetted soil has been considered as the matrix for the bioremediation and phytoremediation processes.OMW degradation in the storage compartment was not considered.
The wastewater's gravitational flow through the soil of the unit was considered to be one-dimensional.Due to the intermittent irrigation of the unit, the wastewater flow was expected to take place both under saturated and unsaturated conditions and was considered to follow a combination of the Darcy-Buckingham model and the "tipping bucket" approach [19,20]; corresponding considerations as well as the wastewater volume balance of each compartment are provided in detail in Appendix A.2.
The total volume of the wastewater in the system was considered to be gradually reduced both through water evaporation and the plants' evapotranspiration; Experimental data used for model calibration and validation were produced during four OMW treatment cycles performed during the summer months in two consecutive years.At the onset of each experimental cycle, both pilots were supplied with diluted OMW in the storage compartment, which subsequently was pumped continuously to the surface of the units for 12 h/d and infiltrated through the soil.Raw OMW had initial average chemical oxygen demand (COD) concentration for Cycles #1-#3 (1st year) and Cycle #4 (2nd year), corresponding to 2700, 6200, 12,400, and 39,300 mg/L, respectively.The dynamics of the most important OMW physicochemical variables, i.e., COD, biochemical oxygen demand (BOD 5 ), total polyphenols (TPh), total nitrogen (TN), and total phosphorus (TP), were investigated.

Model Development of the Phytoremediation Pilot Units
Due to the macroscopic character of the model, several assumptions were made to simplify hydrodynamics and mass transfer phenomena in the pilots.More specifically, the soil layer of the unit was conceptually divided into four equal-volume consecutive compartments (Figure 1b) in a tank-in-series approach [18] to account roughly for the plugflow type of movement of the wastewater liquid phase through the soil.Assumptions on the characteristics of the plants' root system in the two pilots were based on the most relevant available information found in the literature and are described in detail in Appendix A.1.
OMW was provided to the surface of each unit by drip irrigation, where in sandy soil, the vertical water flow was reported to be much faster than the horizontal flow due to the high infiltration rate.Therefore, OMW was considered to be infiltrated only from a specific fraction of each unit's total soil (wetted soil region), shown schematically in Figure 1b.Numerical estimation of the mass of the total wetted soil region was described in detail by Petoussi and Kalogerakis (2022) [4], and its corresponding percent distribution in the unit was considered to be 15%, 25%, 30%, and 30% for the 1st to 4thcompartments, respectively.Wetted soil has been considered as the matrix for the bioremediation and phytoremediation processes.OMW degradation in the storage compartment was not considered.
The wastewater's gravitational flow through the soil of the unit was considered to be one-dimensional.Due to the intermittent irrigation of the unit, the wastewater flow was expected to take place both under saturated and unsaturated conditions and was considered to follow a combination of the Darcy-Buckingham model and the "tipping bucket" approach [19,20]; corresponding considerations as well as the wastewater volume balance of each compartment are provided in detail in Appendix A.2.
The total volume of the wastewater in the system was considered to be gradually reduced both through water evaporation and the plants' evapotranspiration; corresponding considerations are provided in detail in Appendix A.3.Solutes were considered to be distributed between the solid phase (soil) and the liquid phase (soil solution).Moreover, the liquid as well as the soil phase of each compartment was considered to be homogenous.The adsorption of OMW constituents on the soil was modeled by the linear type adsorption isotherm.Convection was considered as the main mass transfer mechanism within the liquid phase due to intense irrigation.To avoid high computational complexity, the transport of solutes within the liquid phase through diffusion was not explicitly modeled.However, in order to account for the effect of the soil water content on the diffusion of soluble OMW substrates in the liquid phase and their availability to microbial biomass, an approximation of the solutes' diffusion in the soil solution was considered, as described in detail in Appendix A.4.The model considerations, particularly those concerning the wastewater oxygenation, are provided in Appendix A.5.

Model State Variables
OMW organic compounds, such as proteins, simple carbohydrates, and lipids, are easily degraded by microorganisms (i.e., labile), while those that are more resistant to decomposition organics, such as lignin, cellulose, and hemicelluloses [21], are partially degraded and transformed slowly [22].Therefore, in the model, OMW total organic matter, determined by COD, was conceptually divided into two organics pools with respect to their biodegradability degree, i.e., readily biodegradable organic matter (RBOM), determined by BOD 5 , and slowly biodegradable organic matter (SBOM), corresponding to the remaining organic fraction (COD\BOD 5 ).Each of the two organic pools was further divided into two pools with respect to the presence/absence of the phenolic character of their compounds.Polyphenolic compounds, a class of OMW organics, are of special interest in OMW management due to their phytotoxic character.As a fact, polyphenolic substances found in OMW have a wide range of molecular weights (MW) [23], from simple monomeric phenols, such as tyrosol and hydroxytyrosol, to high polymerized polyphenols, such as tannins and anthocyanins [24].High MW polyphenols are produced by partial oxidation and polymerization of simple phenols and are generally hardly degradable [25]; thus, their equivalent COD was considered in the model as the phenolic fraction of SBOM, while the equivalent COD of the low MW polyphenols was considered as the phenolic fraction of RBOM.The concentrations of high and low MW polyphenols of OMW (in gallic acid equivalents) were considered as model state variables (C HPh and C LPh ), both comprising the total polyphenol (TPh) pool.
Eventually, OMW organics were considered in the model to be divided in the following pools: (a) non-phenolic SBOM, (b) phenolic SBOM, (c) non-phenolic RBOM, and (d) phenolic RBOM, where pools (a) and (c) corresponded to state variables C S,NPh and C R,NPh , while (b) and (d) were determined as the equivalent COD concentrations of high and low MW polyphenols.
The non-biodegradable fraction of OMW was not considered since OMW is derived from edible matter [26].Note that the SBOM fraction refers to both soluble recalcitrant organics and organics in particulate form.
OMW contains nitrogen mostly in organic form as amino acids, peptides, etc. [27], which, soon after OMW application on the soil, is considered to be mineralized and available for plant nutrition [27].Therefore, OMW total nitrogen (TN) was also considered to comprise two pools, i.e., one of organic and one of inorganic nitrogen compounds, both modeled as state variables (C org−N and C in−N ).
OMW phosphorus is found mainly in the inorganic form of phosphate salts [27][28][29][30], and, thus, only one phosphorus pool was considered as a state variable (C P ).
Additional state variables were considered for the adsorbed mass of non-phenolic SBOM (S ads,S,NPh ), high MW polyphenols (S ads,HPh )), and organic nitrogen (S ads, org−N ) on plants' roots.
Moreover, OMW DO concentration (C O ), the volume of the soil solution in the wetted region of each soil compartment (V k )), and the wastewater volume in the storage compartment were considered as state variables.
The bacterial (X b ) and fungal (X f ) microbial biomass concentration of bulk and rhizosphere soil were also considered as state variables.However, no specialized populations were considered to be involved in specific processes, such as nitrification, ammonification, or phenol degradation.
The general state variables, as described above, were considered for the different compartments of the system (Table 1), and the number of state variables totaled 73.The general mass balance for the soluble constituents (concerning SVs 1 to 8 of Table 1) in the soil solution of the kth unit compartment is given by Equation (1).
where C i,k is the concentration of the ith constituent in the soil solution of the kth compartment (mg/L), V k is the volume of the soil solution of the kth compartment (L), S i,k is the concentration of the ith constituent adsorbed on the soil of the kth compartment (mg/kgsoil), M k is the wetted soil mass of the kth compartment (kg), and R i,k, total = ∑ j (±) R i,k,j is the sum of removal and production rates of the ith constituent in the kth compartment (mg/h) through process j.Rates R i,k,j are presented in Sections 2.2.2 and 2.2.3, and their contributions in the mass balance of each SV is provided in Table A1 of Appendix A.6.For notation simplicity, soil compartment indicator k is omitted in the equations of the processes' rates.
Soil concentration of the ith soluble constituent S i,k , is considered to follow the linear sorption isotherm at every time instant (Equation ( 2)).
where K sd,i (L/kg) is the partition coefficient of the ith constituent.Combining Equation (1) and Equation (2): Particularly for the dissolved oxygen, the closed mass balance of Equation ( 4) is considered only in case of waterlogged conditions, where wastewater oxygen is not replenished from the air in the soil pores; nevertheless, it is consumed as the wastewater infiltrates through the soil.
where j denotes consumption by (a) microbial growth, (b) maintenance, and (c) removal by plants' transpiration.
For the adsorbed organic constituents on plants' root tissue and the microbial biomass, no exchange between soil layers was considered.Therefore SVs X b and X f (g/kg-soil) follow the growth kinetics given by Equation ( 5) and equivalently, the concentration of adsorbed organics on the root tissue of plants where: S ads,i, where i denotes non-phenolic SBOM, high MW phenols, and org-N and follows the kinetics described by Equation (28).

Bioremediation Processes Microbial Growth
Organic compounds of OMW with high energy content induce a rapid short-term increase in the number of soil copiotrophic bacteria or r-selected species which comprise the first colonizers of the newly added organic matter in the soil [31].This phenomenon (priming effect) is reported to last only a period of several weeks after organics application on soil [32], during which time the bacteria assimilate fast RBOM, gaining the advantage over fungi due to their faster growth rate.OMW phenolics, although usually referred as the main factor of OMW antimicrobial properties [33], also include simple molecules of low MW, such as phenolic acids, i.e., tyrosol or hydroxytyrosol, which comprise the majority of the OMW phenolic fraction [34] and have been shown to serve as a carbon source for microbial growth [35][36][37][38][39]. Fungi species are also reported to have an important role in the biodegradation of high-molecular-weight and recalcitrant organic compounds of OMW [5,40,41].
Heterotrophic bacteria and fungi were considered in the model to grow on both the phenolic and the non-phenolic fractions of the labile organic matter (RBOM), following additive kinetics.For the non-phenolic RBOM substrate (sugars, amino acids, etc.), Monod growth kinetics were considered, while for the phenolic fraction of RBOM, Haldane kinetics were followed in order to describe the growth on the inhibitory substrate.Growth limitation by nutrients (nitrogen and phosphorus) was also considered.Bacterial population was considered as facultative aerobic, and fungi were considered as obligatory aerobic, as are most soil fungi species.In aerobic conditions, Monod type oxygen limitation was considered, whereas in anoxic conditions, oxygen limitation followed the inverse Monod model.In anoxic conditions, NO 3 − −N was considered as an additional growth limiting nutrient since NO 3 − −N has the role of an electron acceptor in anoxic bacterial growth and is generally leaching from soil.NO 3 − −N concentration was not considered, however, as a state variable; instead, it was assumed as a constant fraction of inorganic nitrogen concentration.It is noted that the bacterial population is assumed to be capable of growing in both aerobic and anoxic conditions, according to the current oxic state of the soil.This simplification is also made in activated sludge models, ASM1 and ASM2 [42].Decay of the microbial biomass was also taken into account, which led to the replenishment of both the SBOM and organic nitrogen pool.Microbial growth rates are given below in Equations ( 5)- (10): where µ b and d b (1/h) denote the bacterial specific growth and decay rate, correspondingly.Bacterial growth follows additive kinetics, as given in Equations ( 6) and (7).
where µ max,b is the maximum specific bacterial growth rate, and L b corresponds to the bacterial growth limitation factors for the following growth conditions: L b,1 denotes aerobic growth on non-phenolic RBOM; L b,2 denotes aerobic growth on phenolic RBOM; L b,3 denotes anoxic growth on non-phenolic RBOM; and L b,4 denotes anoxic growth on phenolic RBOM.All limitation factors ϕ are shown in Table A2 of Appendix A.6. f NPh and f Ph correspond to the fraction of non-phenolic and phenolic substrate of RBOM pool, respectively, and are considered to weight the contribution of each component of L b to the total rate according to the availability of each substrate.Phenolic RBOM was estimated stoichiometrically as the equivalent BOD 5 of the low MW phenolic concentration expressed in gallic acid equivalents.Moreover, it was implicitly assumed that both aerobic and anoxic growth prevailed at each soil layer according to the current oxygen level of the liquid phase.
The corresponding growth equations for fungi are: where µ f and d f (1/h) denote the specific fungi growth and decay rate, respectively.
where µ max, f is the maximum specific fungi growth rate, and L f are the fungal growth limitation factors, namely, L f ,1 denotes aerobic growth on non-phenolic RBOM and L f ,2 denotes aerobic growth on phenolic RBOM.

Microbial Decomposition and Removal of Organic Compounds
High MW SBOM is assumed to be firstly broken down through extracellular enzymatic decomposition both by bacteria and fungi (Equations ( 11)-( 13)) in order to form bioavailable organic compounds, directly assimilable by microorganisms for growth and maintenance [43].As a part of the total OMW organics, high MW polyphenols are also subject to microbial decomposition and have been reported to be gradually degraded to simple phenolics via first-order enzymatic processes after OMW application on soil [44].
The decomposition of OMW organic matter was herein modeled following the modeling approach of the hydrolysis process in ASM models [42], according to its classical definition, whereby "high molecular weight complex substances are broken down into simpler easily digestible molecules which are subsequently taken up for biomass build up" [42].Enzymatic decomposition was considered, therefore, as surface reaction occurring in close contact between the organisms which provided the hydrolytic enzymes and the slowly biodegradable matter.Moreover, enzymatic decomposition took place in both aerobic and anoxic environments, and although it is not considered to be related to the electron acceptor type, reduction of its rate under anoxic conditions was considered, as it is well described for the anoxic hydrolysis process [42].
Total decomposition rate of the ith compound R i,dec is given as the sum of the corresponding decomposition rates by bacteria R i,b,dec and fungi R i, f ,dec (Equations ( 11)-( 13)), where i accounts for the non-phenolic (i: S, NPh) or phenolic (i: HPh) pool of SBOM.
where k dec, i,b and k dec, i, f (mg-substrate/g-biomass•h) are the maximum enzymatic decomposition specific rates of the ith pool of SBOM for bacteria and fungi, correspondingly, n h (dimensionless) is the reduction factor of enzymatic decomposition rate in anoxic conditions, and are the dimensionless limiting factors for the decomposition process due to the surface limited kinetics.k X corresponds to the saturation constant of the surface limitation terms.
Decomposed non-phenolic SBOM was assumed to contribute equally to the nonphenolic RBOM pool in terms of COD.The same consideration was made for the high MW polyphenols, contributing to the low MW phenolics pool.However, reverse polymerization of simple OMW phenols was also reported to take place [45], which could have possibly led to the underestimation of the real phenolic fraction of the low MW phenolic pool since secondary metabolites of low MW polyphenols can possess higher total antioxidant activity than the initial complex phenolic molecule.However, this could be anticipated by the fact that high MW polyphenols have been also reported to possess high antioxidant activity [46].
Non-phenolic RBOM is consumed during microbial growth according to Equation (14).
where Y X b /R,NPh and Y X f /R,NPh (g-biomass/mg-BOD 5 ) are the bacterial and fungal biomass stoichiometric yields, respectively, on the non-phenolic organic matter.Simple low MW phenolics were also mentioned to be degraded and consumed under aerobic conditions as assimilable carbon substrate from a wide variety of bacteria and fungi species [35][36][37][38][39].However, the inhibition of the polyphenols' biodegradation in pure and mixed cultures at high phenol concentrations is well reported in both aerobic and anaerobic environments [9,47].
Therefore, in the model we considered low MW polyphenol consumption by both bacteria and fungi during microbial growth under inhibition (Equation ( 15)).
where Y X b /LPh and Y X f /LPh (g-biomass/mg-gallic acid) are the biomass yields on low MW polyphenols for bacteria and fungi, correspondingly (Y X y /LPh = Y X y /R,NPh c f ,PhBOD , where c f ,PhBOD is the equivalent BOD 5 of gallic acid).RBOM is also consumed by microbial biomass for maintenance energy (Equation ( 16)), i.e., the energy consumed for functions other than the production of new cell material.
where k m,b and k m, f (mg-substrate/g-biomass•h) are the maintenance coefficients for bacteria and fungi, correspondingly, defined by Pirt (1982) [48] as the minimum substrate consumption to maintain the cells.We also considered simple OMW phenolics to be consumed for microbial maintenance (Equation ( 17)).
A fraction of the decayed microbial biomass, expressed by the dimensionless parameter r d , is considered to contribute back to the non-phenolic SBOM pool (Equation ( 18)).
where c f ,b,SBOM and c f , f ,SBOM are the conversion coefficients of decayed microbial biomass to non-phenolic SBOM based on stoichiometric calculations on the bacteria and fungi chemical formulae, respectively.

Microbial Transformation and Removal of Nitrogenous Compounds
In the model, two nitrogen pools were considered to be contained in OMW, i.e., one pool of inorganic and one of organic nitrogen.The inorganic nitrogen pool was assumed to be enriched with NH 4 + -N by mineralization of organic nitrogen through ammonification (Equation (19).
The ammonification process was taken into consideration in the model, since an increase of the ammonifying bacteria population is often reported in soils amended with OMW [49].Ammonifiers were not modeled, however, as a discrete bacterial population, as was the case in ASM1 model [42].Instead, it was assumed that the possible effect of the different population sizes and kinetics of ammonifiers was lumped in the corresponding process kinetics.Ammonification was modeled as a surface process, similar to organic decomposition (Equation ( 12)): where k amm (mg-substrate/g-biomass•h) is the maximum ammonification specific rate, and A lim,org−N,b (dimensionless) is the surface limitation factor following the same formula as the corresponding factor for SBOM decomposition, as described in the previous section + -N to NO 3 − −N was not modeled, instead NO 3 − −N was assumed to be a constant fraction of the inorganic nitrogen pool at every time instant.Our assumption was based on the results of many studies on OMW direct application on soil, whereby ammonia oxidizing and nitrifying autotrophic bacteria were reported to be slightly inhibited, and nitrification disruption was also reported [50,51].Moreover, OMW was expected to have a low concentration of NO 3 − −N since, at the same time, the stimulation of the soil denitrifier population was reported [50,51], and NO 3 − −N is well known to be leaching in soils.Mineralized nitrogen is considered to contribute equally to the inorganic nitrogen pool and to be thereafter assimilated by bacteria, fungi, and plants.Removal of inorganic nitrogen through incorporation to bacterial and fungal biomass (immobilization) is considered to follow Equation (20).
where Y X b /N and Y X f /N (g-biomass/mg-substrate) are the stoichiometric biomass yields on nitrogen for bacteria and fungi, respectively.Moreover, the organic nitrogen pool is replenished by a fraction (r d ) of the decayed microbial biomass nitrogen (Equation ( 21)).
Denitrification was not considered in the model since its contribution to nitrogen removal is negligible.

Microbial Removal and Precipitation of Phosphorus (P)
Generally, phosphorus can be removed from wastewater either by precipitation and/or adsorption or by microbial uptake [52].OMW phosphorus was considered in the model as consisting of only one pool of available inorganic constituents, which can be assimilated from microbial biomass for build-up (Equation ( 22)).
where Y X b /P and Y X f /P (g-biomass/mg-substrate) are the stoichiometric biomass yields on phosphorus for bacteria and fungi, respectively.Moreover, the chemical precipitation of P in the soil is considered to be a non-reversible first-order kinetics process (Equation ( 23)) [53]: where k prec (1/h) is the specific P precipitation rate.Precipitation, however, is independent of P sorption and desorption to soil, and the latter processes were modeled separately by linear adsorption.Mineralization of organic phosphorus originating from the decomposition of the decayed microbial biomass was not considered to be significant in the short experimental period and, thus, was omitted.

Microbial Consumption of Oxygen
Oxygen is considered to be consumed via aerobic microbial growth (Equation ( 24)) and maintenance (Equation ( 25)): where Y X b /O and Y X f /O (g-biomass/mg-substrate) correspond to the stoichiometric biomass yields on oxygen for bacteria and fungi, respectively.

Phytoremediation Processes Plant Uptake of Organic and Inorganic Compounds
Plants uptake organic contaminants passively from the soil solution during their transpiration [54].The uptake efficiency is affected by the physical and chemical properties of the organic pollutants, i.e., the molecular mass and hydrophobicity, as well as by the biological characteristics of the plants and the environmental media.The molecular mass of organic pollutants has been reported as an important factor influencing the plant-uptake process, whereby compounds with molecular mass below 1000 can be easily absorbed by plant roots [55].Therefore, sugars and proteins of OMW, which have mainly low-molecular mass, and the most abundant low-molecular-mass polyphenols in OMW, i.e., tyrosol and hydroxytyrosol [24], should also be considered to be taken up by plants [46].
Plants have also been known to uptake phosphorus in orthophosphate form [56] and nitrogen in both inorganic and organic form [57]. Inorganic nutrients are known to be taken up actively due to plant-induced gradients that drive them towards the roots [58].
In order to model the uptake of solutes (organics and nutrients) by plants, we herein adopted the approach of transpiration stream concentration factor (TSCF), which was proposed by Briggs et al. (1982) [59] as an indirect measure of the uptake efficiency.It is defined as: TSCF = Concentration in transpiration stream (mg/L)/Concentration of external solution in contact with root tissues (mg/L).According to Manzoni et al. (2011) [58], the plant uptake rate of both organic and inorganic compounds may be approximated simply as proportional to the mean transpiration rate and the compounds' concentration in the soil solution following Equation (26).
where i denotes non-phenolic RBOM, LPh, in-N, org-N, and P. The TSCF has a maximum value of 1 for passive uptake.TSCF factors larger than 1 indicate that nutrients have been moved against their concentration gradient [60].Oxygen has also been considered to be taken up by plants passively, as contained in the evapotranspirated stream of the wastewater.In this type of model, plant uptake is effectively a "sink" process or loss mechanism for chemicals in the soil, and the eventual fate of the compounds within the plant is not investigated.

Accumulation of Recalcitrant Organics to the Plants' Root Zone
High-molecular-mass organic pollutants, which possess strong hydrophobic properties, are considered to be adsorbed on plant roots [55].The potential of a given xenobiotic to accumulate on the plant root is described by the Root Concentration Factor (RCF), defined by Briggs et al. (1982) [59] as the ratio of the organic compound sorbed on the root (mg/kgfresh root tissue) to the compound concentration in the soil solution (mg/L).Therefore, organic accumulation on plants' roots has been modeled according to Equation (27).
where k ads (1/h) is the specific rate of adsorption on roots, and RCF i (mg-substrate/kg-root tissue)/(mg/L) is the root concentration factor of constituent i: non-phenolic SBOM, high MW phenols, and org-N.Rate of change of accumulated mass on the plants' root tissue is given by Equation (28).
Organics while sorbed on roots are not considered to be biodegraded and are subsequently backwashed to the soil solution according to the adsorption potential.

Rhizodegradation of Organic Compounds
Rhizodegradation refers to the breakdown of contaminants within the plant root zone, or rhizosphere.The rhizosphere is a microhabitat around roots, typically extending 1-2 mm around them [61], where an enhanced number and variety of microorganism species (bacteria and fungi) can be found.This "rhizosphere effect" is caused by the physical (i.e., gas exchange and soil moisture) and chemical (root exudates) impact of plant roots on the soil [54].Rhizosphere microbial population concentration generally varies from 10-to 1000-fold with respect to bulk soil microbial biomass concentration [62].Although no experimental assessment of the rhizosphere microbial population in the pilots was available, rhizodegradation could not be neglected in the model.Therefore, the rhizosphere microbial biomass concentration was modeled assuming additional state variables to describe the bacterial and fungal population in each of the 4 different soil compartments.The initial rhizosphere fungal and bacterial biomass concentration in the pilots was assumed to be two orders of magnitude greater than the corresponding concentration of the bulk soil.Since this assumption was highly uncertain, it was further examined in the sensitivity analysis of the model.It is noted also that rhizosphere bacterial and fungal biomass were considered to follow the same bioprocesses and relevant kinetics with the corresponding bulk populations.The mass of the total rhizosphere soil in the experimental period of the first year was estimated as 10% of the total wetted soil region of each pilot (20% for second year); this was estimated according to the calculation of the soil volume surrounding roots based on available information about the plants' root length distribution [63]) and the corresponding rhizosphere activity length [64].The total mass of rhizosphere soil was distributed to the four unit layers according to the root distribution profile of each plant species.No additional rhizosphere soil compartments were considered in the model, for simplicity reasons.Due to the fast wastewater infiltration through the soil, it was assumed that the concentration of solutes in the soil solution of the bulk and the rhizosphere soil were the same.

Model Inputs and Parameters
The parameters of the model to be calibrated are related to the processes of (a) plant uptake, (b) adsorption on plants' roots, (c) microbial growth, (d) enzymatic decomposition of organics, and (e) phosphorus precipitation.These are correspondingly (a) transpiration stream concentration factors TSCF i (i denotes non-phenolic RBOM, low MW phenolics, in-N, org-N, and P), (b) root concentration factors RCF i (i denotes non-phenolic SBOM, high MW phenolics, and org-N) and adsorption-specific rate constant k ads , (c) maximum specific microbial growth rates µ max,y (y denotes bacteria and fungi), semi-saturation constants K S,i,y, (i denotes non-phenolic RBOM, phenolic RBOM, in-N, and P; and y denotes bacteria and fungi), semi-saturation constant K S,NO3,b for NO 3 − N in anoxic growth conditions, and inhibition constants for phenolic RBOM K S i ,R,Ph,y (y denotes bacteria and fungi), (d) factor of enzymatic decomposition rate reduction in anoxic conditions n h , enzymatic decomposition specific rates k dec,i,y (i denotes non-phenolic SBOM and high MW phenolics; and y denotes bacteria and fungi), and ammonification specific rate k ammon , and (e) phosphorus precipitation specific rate k prec .Units of the model inputs and the kinetic parameters to be calibrated are provided, correspondingly, in Tables A3 and A4 of Appendix A.6.

Model Calibration Methodology
The mechanistic model of the present study was developed in MATLAB ® in statespace form.The governing differential equations were solved by MATLAB's ode15s differential equation solver.The model calibration was based on the available experimental data series of BOD 5 , COD, TPh, TN, and TP concentrations of the recirculating OMW from Treatment Cycles #1 and #4 (i.e., treatment of low-and high-organic-strength OMWs, correspondingly) as well as the data of soil concentrations of bacterial biomass, total organic carbon (TOC), and TN measured at the end of the second experimental year given by Petoussi and Kalogerakis (2022) [4].The corresponding data from Cycles #2 and #3 (treatment of medium-organic-strength OMW) were used for model validation.Due to the different soil environments created by each plant species, the calibration procedure was conducted separately for the two phytoremediation pilot units.The calibration was performed according to the adaptive LJ optimization procedure [65], which uses random search points and systematic contraction of the search region.According to this method, an initial guess k(0) of the p-dimensional unknown parameter vector (p = 29) was made and the initial search region was set as r(0) = k max − k min .The sum of the normalized root mean square errors (NRMSEs) of the experimentally defined variables was considered as the objective function to reflect the agreement of model prediction with measured data (Equation ( 29)): where N pm is the total number of experimentally measured variables, N j is the total number of experimental data points available for the jth variable, O ij and P ij are the measured and corresponding predicted values, respectively, by the model for the ith experimental data point of the jth variable, and O j is the mean value of the measured experimental data points.The number of random evaluations of the objective function within an iteration was set to 100, while the maximum number of iterations was set to 40, with a contraction of the search region at the end of each iteration by 95%.

Simulated Dynamics
The simulated dynamics of OMW organic concentration during Cycle #1 and Cycle #4 for the two experimental units are shown in Figure 2.
The model follows satisfactorily the datasets of COD (Figure 2C,G) and TPh (Figure 2D,H) concentrations for both units during Cycle #4 (NRMSE equals 24% and 8% for P. granatum and M. communis pilots, correspondingly), where OMW of high organic strength is treated.On the contrary, during Cycle #1, COD and TPh concentrations are predicted with a high bias (average NRMSE for both units is equal to 78%) from the corresponding datasets (Figure 2A,B,E,F).The divergence of the model prediction from the data is due primarily to the limitation of the microbial growth by the lack of nutrients since labile organic matter is available throughout the entire experimental period of Cycle #1.During the last days of Cycle #1, the model also tends to increase the COD concentration on the basis of (a) the assumption of redissolution of organic matter concentrated in the plants' roots back to the soil solution and (b) the replenishment of the SBOM pool by the decayed biomass.As a fact, the model predicts less accurately the system dynamics during the treatment of low organic strength wastewater.From a process-related view, this also could be attributed to the promotion of organics enzymatic degradation in highly diluted OMW [25], which is not included mechanistically in the model.
Data regarding TN and TP concentrations are generally in high agreement with the model dynamics (Figures S3 and S4) for both units and cycles.

Simulated Dynamics
The simulated dynamics of OMW organic concentration during Cycle #1 and Cycle #4 for the two experimental units are shown in Figure 2. The simulated dynamics of soil bacterial biomass, TOC, and TN concentrations are given in Figure 3 for both the bulk and rhizosphere soil for M. communis (Figure 3A-C) and P. granatum units (Figure 3D-F) during Cycle #4.
Similar dynamics for rhizosphere and bulk soil microbial populations result from the fact that the same substrate availability and growth kinetics have been considered for the two populations.Bacterial biomass concentration (Figure 3A,D) increases initially due to the soil priming effect caused by the application of OMW organic matter, as also described by Kotsou et al. (2004) [31].However, although labile organic matter is highly available up to the end of Cycle #4 (Figure 2C,G), the bacterial biomass concentration, especially in the P. granatum unit soon after the beginning of the cycle, follows a decreasing trend.This is due to the simultaneous substrate consumption by fungi and to the multiple nutrient limitation of the microbial growth.At the end of Cycle #4, there is a perfect match of the model prediction with data for the bacterial biomass in both units.Simulated dynamics of fungal biomass concentration are given as Supplementary Materials (Figures S8 and S9) since no corresponding experimental data were available for model calibration or validation.
[25], which is not included mechanistically in the model.
Data regarding TN and TP concentrations are generally in high agreement with the model dynamics (Figures S3 and S4) for both units and cycles.
The simulated dynamics of soil bacterial biomass, TOC, and TN concentrations are given in Figure 3 for both the bulk and rhizosphere soil for M. communis (Figure 3A-C) and P. granatum units (Figure 3D-F) during Cycle #4.Similar dynamics for rhizosphere and bulk soil microbial populations result from the fact that the same substrate availability and growth kinetics have been considered for the two populations.Bacterial biomass concentration (Figure 3A,D) increases initially due to the soil priming effect caused by the application of OMW organic matter, as also described by Soil TOC and TN concentrations (Figure 3B,C,E,F) acquire and maintain higher values in the rhizosphere soil (which computationally include the accumulated recalcitrant organics on plants' roots) as compared with the bulk soil in both pilots.Specifically in the case of the M. communis plants pilot, the rhizosphere soil acquires disproportionally high TOC and TN concentrations (Figure 3B,C), likely due to the root distribution considered for M. communis plants, where the highest proportion of the roots' biomass (79%) is considered to be developed in the top soil compartment.
Experimentally defined values of TOC and TN concentrations for the four soil compartments are not significantly different in any pilot unit, as determined by Petoussi and Kalogerakis (2022) [4].In addition, there is no certainty that the soil samples were taken from the bulk soil excluding rhizosphere soil.Therefore, the model prediction is satisfactory even in the M. communis case since data points of TOC and TN concentrations at the end of the cycle lie between the predicted concentrations of rhizosphere and bulk soil (Figure 3B,C).

Parameter Estimation
Using the LJ optimization procedure, we arrived at the estimated parameter values for the two experimental units shown in Table A4 of Appendix A.6.

Kinetics of Microbial Growth
Maximum specific microbial growth rates for bacteria and fungi (µ max,b and µ max, f ) were estimated by the model in the narrow range of 0.164-0.2361/h for both pilots.The exception of µ max,b for the M. communis unit estimated as 0.027 1/h (one order of magnitude lower) was examined in the sensitivity analysis.Corresponding specific growth rates reported in studies of OMW degradation kinetics in bioreactors are generally up to two orders of magnitude lower both for bacteria [9,15,17,66] and fungi [5,15,40,67,68], irrespective of the kinetics model used.This inconsistency could result from the fact that in the model developed in this study, multiple substrate limitation was considered, which required higher maximum growth rates.
Semi-saturation constants for the microbial growth on non-phenolic substrate K S,R,NPh,b and K S,R,NPh, f take values in the range of 32.26-101.19mg/L, which is lower than the corresponding range of 232-860 mg/L reported by Chiavola et al. (2014) [9], who modeled the kinetics of raw OMW aerobic treatment by activated sludge biomass according to the Monod model.Much higher semi-saturation constants in the range of 4000-18,000 mg/L have also been reported [40,68] for the Monod kinetics growth of fungi on raw OMW.For the microbial growth on low MW phenolics, the semi-saturation constants K S,R,Ph,b and K S,R,Ph, f are in the range of 56.63-82.01mg/L,and inhibition constants K S i ,R,Ph,b and K S i ,R,Ph, f are in the range of 23.12-129.87mg/L.Similar values of the previous constants were also reported by Tziotzios et al. (2008) [66],who modeled phenol removal by olive pulp bacteria growing in a batch culture following Haldane growth kinetics.The inhibition constants estimated in our model are remarkably low, which gives evidence of the inhibitory action of the phenolic organic substrate [66].

Kinetics of Organic Enzymatic Degradation
Specific rates k dec, S,NPh for the decomposition of SBOM are in the range of 10.75-72.73mgsubstrate/g-biomass•h both for bacteria and fungi.Corresponding rates for high MW phenols decomposition k dec,HPh are in the range of 1.02-18.25 mg-substrate/g-biomass•h.The most relevant literature for the decomposition of organic matter in the soil concerns studies in composting of organic wastes, where organic decomposition has been generally modeled as a first-order kinetics process.Corresponding kinetic constants reported for the composting of OMW [16], urban organic waste [69], and specialized organic pollutants, such as PAH, surfactants, and herbicides [70], are more than three orders of magnitude smaller than the estimated kinetic values of our model.The high kinetic constants of this study likely counterbalance the surface limitation considered for the enzymatic decomposition.

Kinetics of Phytoremediation Processes
The estimated value of the transpiration stream concentration factor for the nonphenolic organics TSCF R,NPh equals 1 for both units, resulting from the imposed constraint on the maximum theoretical value of TSCF for the passive plant uptake along with the water used for transpiration.Experimentally defined values of TSCF greater than unity have already been reported in the literature [71]; however, they are not considered to be data of high significance [72].Values of TSCF for low MW phenols are also greater than 0.95.It was expected that TSCF values would be different for the two plants since uptake efficiency has been reported to vary with plant species, age, and health among others [72,73].TSCF values derived from the model calibration cannot be directly compared with the corresponding values of the various chemicals reported in the literature [71] since OMW is a mixture of many substances, and TSCF in this case refers to a class of constituents and not to one chemical substance.In addition, large variations have been observed in the literature for compounds with more than one average TSCF, which was attributed to the different experimental approaches and operational variables [74].TSCF values of inorganic fractions of in-N and P are several times greater than unity, with the exception of the very low phosphorus TSCF value (0.02) for P. granatum plants (examined in the sensitivity analysis in Section 3.4).The TSCF value can be greater than unity for nutrients such as nitrogen and phosphorus if they are actively taken up by the plants [75].Values of TSCF for inorganic nutrients are not available in the literature.
Root concentration factor (RCF) values for the recalcitrant OMW organics (nonphenolic SBOM, high MW phenols, and org-N) in the two units are within the range of the reported values in the relevant literature [76].Namiki et al. (2018) [77] discussed the different abilities among plant species to take up organic chemicals from the soil into their roots, and they reported RCFs of organic chemicals taking a high range of values.

Model Validation
As shown in Figure 4, the model predicts very well the dataset of COD (Figure 4C,G) and TPh concentrations (Figure 4D,H) during Cycle#3 for both pilot units (average NRMSE equals 27% and 31% for COD and TPh, correspondingly).Soil TOC, TN, and bacterial biomass concentrations, as shown in Figure 5 for the entire experimental period of the 1st year (Cycles #1-#3),are increasing at the onset of each OMW treatment cycle, as fresh OMW is applied to the soil (time instants t=0 d, t=15 d, and t=27 d).During Cycle #2, there is an adequate agreement of the COD (Figure 4A,E) as well as the TPh (Figure 4B,F)concentration data with the model prediction, albeit only until the first half of the cycle.After the 6th day, a positive bias in the prediction of the model is present in both pilots.TN (Figure S5) and TP (Figure S6) concentration data are, however, predicted accurately for Cycles #2 and #3 for both plant units at all time instances, with the exception of TP prediction during Cycle #3, where the data present high fluctuation which cannot be predicted by the model (Figure S6b,d).
Soil TOC, TN, and bacterial biomass concentrations, as shown in Figure 5 for the entire experimental period of the 1st year (Cycles #1-#3),are increasing at the onset of each OMW treatment cycle, as fresh OMW is applied to the soil (time instants t = 0 d, t = 15 d, and t = 27 d).Soil TOC, TN, and bacterial biomass concentrations, as shown in Figure 5 for the entire experimental period of the 1st year (Cycles #1-#3),are increasing at the onset of each OMW treatment cycle, as fresh OMW is applied to the soil (time instants t=0 d, t=15 d, and t=27 d).The bacterial growth rate is proportional to the organic load of the wastewater, as is also referenced by Magdich et al., (2013) [78] and Mekki et al. (2018) [79].The model prediction is satisfactory since data points lie between the predicted concentration of rhizosphere and bulk soil, with the exception of the TN prediction at the end of the experimental year in the M. communis unit, which has a relatively low NRMSE of 21%.

Contribution of the Phytoremediation Processes in the Removal of OMW Organics
Figure 6 shows the relative contribution of the various removal pathways for organics and TPh in the M. communis and P. granatum pilot units.The bacterial growth rate is proportional to the organic load of the wastewater, as is also referenced by Magdich et al., (2013) [78] and Mekki et al. (2018) [79].The model prediction is satisfactory since data points lie between the predicted concentration of rhizosphere and bulk soil, with the exception of the TN prediction at the end of the experimental year in the M. communis unit, which has a relatively low NRMSE of 21%.

Contribution of the Phytoremediation Processes in the Removal of OMW Organics
Figure 6 shows the relative contribution of the various removal pathways for organics and TPh in the M. communis and P. granatum pilot units.
Organic removal pathways considered in the model are (a) removal through plant uptake and (b) removal through microbial growth and maintenance in rhizosphere and bulk soil.Mass not removed from the system refers to the total organic mass distributed in the system, e.g., sorbed on the soil, dissolved in the soil solution, concentrated on the plants' roots, or dissolved in the stored wastewater.Corresponding results for TN and TP removal are provided as Supplementary Materials (Figure S7).
As shown in Figure 6, the main organic removal pathway for both units during Cycle #1 is plant evapotranspiration (26-33% COD removal and 18-22% TPh removal).However, during the next cycles (#2 and #3), plant evapotranspiration contributes less to organic removal, whereas consumption by rhizosphere microbial biomass gradually increases its contribution up to 5-8% COD and 4-10% TPh removal at the end of Cycle #3, where more bioavailable organic substrate and nutrients are present in the wastewater, thus promoting microbial growth.At Cycle #3, plant evapotranspiration and consumption by rhizosphere microbial biomass contribute almost evenly to organic removal and together by 11-13% COD removal and 12-18% TPh removal.Plant evapotranspiration is the main pathway for organic removal at Cycle #4, corresponding to 25-29% COD and 12-15% TPh removal, which could be attributed to the higher evapotranspiration rate of plants during the 2nd experimental year.
diction is satisfactory since data points lie between the predicted concentration of rhizosphere and bulk soil, with the exception of the TN prediction at the end of the experimental year in the M. communis unit, which has a relatively low NRMSE of 21%.

Contribution of the Phytoremediation Processes in the Removal of OMW Organics
Figure 6 shows the relative contribution of the various removal pathways for organics and TPh in the M. communis and P. granatum pilot units.Organic removal pathways considered in the model are (a) removal through plant uptake and (b) removal through microbial growth and maintenance in rhizosphere and bulk soil.Mass not removed from the system refers to the total organic mass distributed in the system, e.g., sorbed on the soil, dissolved in the soil solution, concentrated on the plants' roots, or dissolved in the stored wastewater.Corresponding results for TN and TP removal are provided as Supplementary Materials (Figure S7).
As shown in Figure 6, the main organic removal pathway for both units during Cycle #1 is plant evapotranspiration (26-33% COD removal and 18-22% TPh removal).However, during the next cycles (#2 and #3), plant evapotranspiration contributes less to or- TN (Figure S7a,b) is removed mainly through plants uptake and by more than 79% and 61% in M. communis and P. granatum pilots, respectively, in all cycles of the first year.Computationally, this performance is in line with the active uptake of inorganic nitrogen by plants and in accordance with the corresponding TSCF values greater than unity (Table A4).However, at Cycle #4, consumption by rhizosphere microbial biomass is also contributing highly to TN removal by 50% and 20% of the total wastewater TN mass for the P. granatum and M. communis units, respectively.It is noted that denitrification was neglected in the model for the purpose of simplification, as it contributed less than 1% in the removal of TN.
In the case of TP (Figure S7c,d), there is a difference between the two pilots regarding the main removal TP pathway, i.e., in the P. granatum unit, phosphorus is precipitated mainly in all cycles (>86% in 1st year's cycles and 61% in Cycle #4), whereas in the M. communis unit, it is also taken up by plants in a high percentage (27-40% in 1st year's cycles and 45% in Cycle #4 of the 2nd year).This is computationally attributed to the large discrepancy between the TSCF values for TP in the two units (TSCF P equal to 0.02 and 2.83 for the P. granatum and M. communis units, correspondingly) since precipitation constant values k prec are similar for the two pilots.
The contribution of the bulk soil microbial biomass activity to the removal of all wastewater compounds is minimal (less than 1%) in all cycles for both units, and more specifically, it does not exceed 9% of the total removed organic amount in any of the cases.Overall, this indicates that the phytoremediation processes (plant uptake and rhizodegradation) contribute in all cases by more than 91% to the organic removal in the context of the present study.It is noted, however, that removal percentages should not be compared among experimental cycles, as they refer to different initial wastewater concentrations, i.e., the decreasing trend of the removal of organics from Cycle #1 to Cycle #3 is not indicative of a less effective system.
It is also noted that the mass of the removed organics (COD and TPh) from the pilot system estimated by the model is, in all cases, substantially lower than the one evaluated by Petoussi and Kalogerakis (2022) [4] since the model considers accumulation on the plants' root tissue in addition to adsorption on soil.

Sensitivity Analysis
A sensitivity analysis of the model was performed according to the procedure outlined by Zhang et al. (2012) [69].Each model parameter was increased independently by 5% from its estimated value.The sum of the deviation of all SVs considered for the wastewater in the storage compartment and SVs of microbial biomass concentration at specific time instances of each of the four experimental cycles were considered as the sensitivity coefficient σ (Equation ( 30)).
SV(P 0 + 5%, t) − SV(P 0 , t) SV(P 0 ) t = 0, 1, 2, . . ., N f (30) where N cycles is the number of experimental cycles (4), N sv is the number of SVs considered for the sensitivity analysis, SV(P 0 , t) is the value of the state variable SV at time t when parameter P has its calibrated value P 0 , SV(P 0 + 5%, t) is the value of the state variable SV at time instant t when parameter P has its calibrated value Po increased by 5%, and SV(P 0 ) is the mean value of the SV.The sensitivity coefficients of all estimated model parameters are presented in Table A5 of Appendix A.6 for both pilot units.Parameters with sensitivity coefficient σ > 0.01 were further assessed for their impact on the system's removal efficiency for COD and TPh based on the simulation of Cycle #4 (with high-organic-strength wastewater).In this context, each parameter was increased by a certain percentage, and the normalized percent variation in the removal of COD and TPh by each process was assessed.A few additional model inputs related to the plants and wastewater characteristics were also investigated.Results are provided in Table A6 of Appendix A.6.As shown, the estimated parameters causing the greatest percent variation in COD and TPh removal (with absolute coefficients greater than 0.3) are related to nutrients availability and maximum microbial biomass growth specific rates and, therefore, have an impact on organic removal through microbial consumption (biomass growth and maintenance).Among those parameters, k prec and µ max, f have the greatest impact on the removal of COD and TPh for the P. granatum and M. communis units, respectively.However, for a 20% increase in those parameters, an increase in COD and TPh removal does not exceed 8% and 25%, correspondingly, for the two pilot units, which, in both cases, provides no more than 1.5% additional removal.
Among the parameters related to the wastewater synthesis, the percentage of high MW polyphenols to total polyphenols C HPh C HPh +C LPh 0 affects strongly affects TPh removal through both microbial consumption and plant evapotranspiration in both units.This is expected, due to the change of the low MW phenols availability for the previous processes.A 20% increase in high MW polyphenols leads to a reduction in TPh removal up to 57% through plant uptake and up to 47% through rhizosphere microbial biomass activity in both units.
In addition, important variation in organic removal is also caused by the parameter m tissue related to the root system of the plants, which is expected since it causes variation in the concentration of organics on the plants' roots.A 20% increase in m tissue parameter leads to a 9% and 14.5% reduction in COD and TPh removal, correspondingly.
Concerning the kinetic parameters which are estimated with very diverse values for the two pilot units (Table A4), the sensitivity analysis shows that µ max,b , TSCF P , and k S,in−N, f are more sensitive in the M. communis pilot system, and k S,P, f and K S,in−N,b are more sensitive in the P. granatum system; thus, the estimated values for the corresponding systems are considered to be more reliable.
Generally, according to the sensitivity analysis results, the removal efficiency of the system is expected to be improved in the case of an extra supply of nutrients N and P in a soil environment where the growth of fungi species is promoted.For higher efficiency of the system in TPh removal, treatment of the wastewater soon after its production should be considered in order to maintain a lower fraction of high MW polyphenols.
The developed model is macroscopic and is based on available experimental data [4].Therefore, it has some limitations, among which the most important are: (a) the effect of OMW pH and temperature on the rates of the microbial processes was taken into account, (b) the modeling of highly computationally demanding diffusion processes, as well as the consideration of additional soil compartments for the rhizosphere zone, were excluded, which could affect the availability of OMW solutes for both microbial consumption and removal by plants through transpiration, (c) the model calibration was not based on the soil fungal concentration (due to lack of experimental data), which could underestimate the role of the corresponding microbial population in OMW treatment, (d) the accumulated recalcitrant organics on the plants' roots were not considered to be degraded in the interest of model simplification, which is not very realistic since they are degraded by the rhizosphere microbial biomass, and (e) the values of the roots' tissue mass of the plants m tissue as well as the OMW initial high MW polyphenols fraction

Conclusions
An integrated macroscopic mechanistic model was developed to describe the dynamics of OMW degradation during its short-term recirculation in the soil of two phytoremediation pilot units using plants of the P. granatum L. and M. communis L. species.An innovation of the model is the discretization of the total labile organic substrate of OMW to a phenolic and a non-phenolic fraction, consumed by both bacteria and fungi.The model effectively predicted the experimental data from the treatment of medium-and high-organic-strength wastewater.The maximum specific kinetic rates of microbial growth and organics decomposition were evaluated through the model calibration process to be more than two orders of magnitude higher than the corresponding rates reported in studies of OMW degradation in bioreactors or in studies of composting processes for various organics.The values of TSCF for the non-phenolic and phenolic fractions of OMW organics were estimated as greater than 0.95, which indicate almost passive uptake of the OMW organic mixture by the plants.TSCF factors of inorganic N and P have been estimated as greater than unity, indicating corresponding active uptake of these ingredients.
The relative contribution of plant uptake and rhizodegradation processes in the OMW organic removal varies among the treatment cycles since it depends on both the evapotranspiration rate of plants and the microbial biomass concentration in their rhizosphere.However, the total contribution of both processes on the overall treatment of OMW was specified as more than 91% during treatment of low-or high-organic-strength OMW in both pilot units.
The estimated kinetic parameters causing the greatest percent variation in organic removal are related to the nutrients availability in the system and the microbial maximum specific growth rates, as investigated through the sensitivity analysis of the model.However, for a 20% increase in those parameters, the additional absolute removal of OMW organics mass did not exceed 1.5%.The most important parameter for the removal of polyphenols was the percentage of high MW polyphenols concentration over the total polyphenols concentration.An increase in the previous parameter caused a substantial decrease in polyphenols removal from the system.
Generally, according to the sensitivity analysis results, the efficiency of the system is expected to be improved in the case of an extra supply of nutrients N and P to the system and promotion of the fungal growth.For higher efficiency of the system in TPh removal, treatment of the wastewater soon after its production should be considered in order to maintain a lower fraction of high MW polyphenols.Symbol "*" denotes the states for the solute concentration at the biofilm interface C int (Appendix A.4). k S,i,y (mg/L) denotes semi-saturation constant of constituent i in the growth rate of population y. i denotes the states for non-phenolic RBOM, phenolic RBOM, in-N, P, and oxygen; and y denotes the states for bacteria or fungi; K S i ,R,Ph,y (mg/L) denotes inhibition constant for phenolic RBOM in the growth rate of population y; η NO 3 (dimensionless) denotes NO 3 − N fraction of in-N.

Figure 1 .
Figure 1.(a) Schematic view of the phytoremediation pilot units and (b) compartments considered by the model.

Figure 1 .
Figure 1.(a) Schematic view of the phytoremediation pilot units and (b) compartments considered by the model.

Figure 2 .
Figure 2. Dynamics of COD, BOD 5 , and TPh concentrations of OMW in the storage compartment of (A-D) P. granatum and (E-H) M. communis plants unit during Cycles #1 and #4 of the experiment.

Figure 3 .
Figure 3. Dynamics of the bacterial biomass, TOC, and TN concentrations of the soil in (A-C) P. granatum and (D-F) M. communis plants units during Cycle #4 (2nd year) of the experiment.

Figure 3 .
Figure 3. Dynamics of the bacterial biomass, TOC, and TN concentrations of the soil in (A-C) P. granatum and (D-F) M. communis plants units during Cycle #4 (2nd year) of the experiment.

Figure 4 .
Figure 4. Dynamics of COD, BOD5,and TPh concentrations of OMW in the storage compartment of (A-D) P. granatum and (E-H) M. communis plants units during Cycles #2 and #3 of the experiment.

Figure 4 .
Figure 4. Dynamics of COD, BOD 5 ,and TPh concentrations of OMW in the storage compartment of (A-D) P. granatum and (E-H) M. communis plants units during Cycles #2 and #3 of the experiment.

Figure 4 .
Figure 4. Dynamics of COD, BOD5,and TPh concentrations of OMW in the storage compartment of (A-D) P. granatum and (E-H) M. communis plants units during Cycles #2 and #3 of the experiment.

Figure 5 .
Figure 5. Dynamics of the bacterial biomass, TOC, and TN concentrations of the soil in (A-C) P. granatum and (D-F) M. communis plants units during Cycles #1-3 (1st year) of the experiment.

Figure 6 .
Figure 6.Percentage contribution of the various processes for (A,B) COD and (C,D) TPh mass removal in P. granatum and M. communis units for all experimental cycles.

CDFigure 6 .
Figure 6.Percentage contribution of the various processes for (A,B) COD and (C,D) TPh mass removal in P. granatum and M. communis units for all experimental cycles.
LPh 0 were based on the relevant literature data (due to lack of experimental data).

Table 1 .
List of state variables used in the model.

Table A1 .
Kinetic rates of all processes in the mass balance of each state variable.

Table A2 .
List of limitation expressions for microbial growth.

Table A4 .
Values of the calibrated parameters of the model.(mg-substrate/kg-root tissue)/(mg/L), ** mg-substrate/g-biomass h, *** b and f indicators in parameter symbols account for bacterial and fungal biomass, correspondingly. *

Table A5 .
Sensitivity coefficients of σ of the calibrated parameters for both pilot units.