Green Roofs for domestic wastewater treatment: Experimental and numerical analysis of nitrogen turnover

The use of Green Roofs (GRs) for domestic wastewater treatment and nitrogen (N) removal is an appealing opportunity to conjugate hydrological, energetic, and water quality benefits. However, the research in this direction has been limited to few experimental studies, while the role of numerical modeling for analysis and design has been overlooked. To advance understating of reactive transport processes in GRs and investigate the use of numerical models, this study presents a comprehensive experimental and numerical investigation of the N turnover in vegetated and non-vegetated GRs exposed to domestic wastewater injection. A Bayesian framework is used to calibrate a mechanistic numerical model against multiple observations from a controlled laboratory experiment. A Global Sensitivity Analysis complements the framework to further enlighten the role of multiple physical processes. Results indicate that GRs have high nitrification capacity, which increases after wastewater application. The total leached N is 94% of the total injected N in a non-vegetated GR, compared to 67% for a vegetated GR. Nitrate leaching is negatively correlated with the root solute uptake capacity of the vegetation, which may reach up to 32% of the total injected N. Maximization of this factor is desirable and can be achieved by selecting appropriate plant species and reducing water stress periods using a proper irrigation schedule.


Introduction
In a world with a growing population with increasing demand for food and water, current and future generations face multiple social, economic, and environmental challenges. According to the United Nations (2018), the urban population is projected to increase by 2.5 billion between 2018 and 2050, reaching approximately 68% of the total world population. If not managed properly, this growth would increase the environmental footprint of urban areas by aggravating the already existing hydrological, energetic, and water quality issues. As reported by the European Environmental Agency (EEA), between 20% and 25% of rivers and coastal water bodies are subject to point source pollution from industrial facilities, sewage systems, and wastewater treatment plants (WWTPs) EEA, 2018.
Excess nitrogen input into ecosystems has an extensive detrimental effect on the environment and human health (Hill, 2010;Holmes et al., 2019;World Health Organization, 2011;Yamashita & Yamamoto-Ikemoto, 2014). When untreated, elevated concentrations of ammonia, nitrite, and nitrate in water bodies contribute to eutrophication, leading to hypoxic "dead zones". Unfortunately, many conventional WWTPs are not able to completely remove nitrogen species from wastewater prior to their disposal into water bodies (Holmes et al., 2019;Johnston et al., 2019). Although alternative technologies are highly efficient (e.g., membrane technologies), they are usually expensive and energy-intensive, thus hindering their widespread adoption (Matassa et al., 2015;Mehta et al., 2015). Consequently, urban rivers have become pools of reactive nitrogen and hotspots of regional pollution .
Green infrastructure represents a potential alternative solution to manage excess nitrogen at the local scale through wastewater recycling and treatment (Yu et al., 2019). The combination of physical, chemical, and biological processes occurring in such infrastructure systems, such as soil water retention, sorption, phytoremediation, and microbial degradation (Dai et al., 2017;Papaevangelou et al., 2016;Wirasnita et al., 2018), can be used to provide a preliminary wastewater treatment diffusively in the urban area (Masi et al., 2017). Considering that rooftops may represent up to 50% of impervious surfaces in urban areas and that domestic wastewater is an important source of nitrogen pollution, the use of Green Roofs (GRs) for wastewater treatment represents an appealing but unexploited opportunity to conjugate multiple environmental benefits (Pradhan et al., 2019). In particular, such systems could be able to 1) retain, delay, and evapotranspirate stormwater (Baek et al., 2020;Berretta et al., 2014;Brunetti et al., 2016a;Li and Babcock, 2014;Stovin et al., 2012) 2) reduce the energy consumption of the building by maximizing the evaporative cooling effect induced by the continuous feeding of domestic wastewater on the GR (Brunetti et al., 2018;Castleton et al., 2010;Lazzarin et al., 2005), and 3) diffusively treat wastewater through a combination of sorption, microbial degradation, and plant's uptake and detoxification processes.
Despite many potential benefits, research on this topic has been limited (Vo et al., 2019). Only a few studies have investigated the performance of GRs for wastewater treatment, usually with promising outcomes. For example, results reported in Zapater-Pereyra et al. (2016), Zapater-Pereyra et al. (2013) indicated high ammonium, total nitrogen, and total phosphorous removal efficiencies of a 9-cm deep GR, mainly due to its high aeration capacity. Similar findings were obtained by Van et al. (2014), Song et al. (2013), and Thanh et al. (2014). Furthermore, there is a lack of studies dealing with the development and testing of numerical models for such systems' analysis and designs. Advancement in the understanding of the physicochemical processes in GRs would clarify their potential benefits and allow for the development of numerical models for design and analysis.
Therefore, the main objectives of this study are to experimentally quantify the nitrogen removal efficiency of GRs treated with wastewater, and numerically describe and investigate the factors driving the nitrogen reactive transport in GRs. To this aim, an experimental campaign involving the application of wastewater on both vegetated and non-vegetated GRs under controlled laboratory conditions was carried out. To characterize water flow and reactive nitrogen transport in the GRs, the experiment was complemented with tracer applications, subsurface outflow and soil moisture monitoring, and nitrogen concentration measurements. Observations were combined with a Bayesian inference framework to calibrate the hydrological model HYDRUS (Šimůnek et al., 2016) and inversely estimate soil hydraulic and solute transport parameters. Finally, the calibrated model is coupled with a Global Sensitivity analysis to enlighten the role of different physical processes on nitrogen turnover and nutrients leaching from the GR.

Experimental setup
The experimental setup ( Fig. 1) included three 15 cm-deep GRs testbeds assembled in plastic boxes (56×37 cm) and equipped with an irrigation system consisting of multiple drippers. One box served as the control, and its hydraulic behavior was already investigated in Brunetti et al. (2020a). The boxes were placed on a wooden frame to insulate the bottom and minimize heat transfer from the surrounding environment. A slope of 1% was set to convey subsurface outflow to a circular hole and then to a plastic pipe for measurement and sampling. The GRs' morphology and the irrigation system were identical in all boxes. A 12cm deep granular substrate, consisting of a mixture of volcanic minerals with a measured porosity and bulk density of 0.66 and 0.9 g/cm 3 , respectively, pH equal to 7, initial nitrogen and organic matter contents of 4 and 8% (Brunetti et al., 2016a), respectively, was used as a growing medium. A highly permeable geotextile separated the soil substrate from the underlining 3 cm-deep drainage layer composed of coarse gravel. Two testbeds were vegetated, and one was left with a bare soil surface. Different designs were meant to highlight the role of vegetation on the coupled water-solute transport. Two plant species frequently used for extensive green roofs were installed (Brunetti et al., 2016a): Lobularia Maritima and Dianthus Gratianapolitanus. Both plants were transplanted after a two-week equilibration period, during which only water was applied on the testbeds using peristaltic pumps.
Testbeds were equipped with two thermocouples to measure soil temperature at two different depths and a coaxial impedance dielectric reflectometer (Stevens HydraProbe) at z = − 6 cm to measure soil moisture and characterize the hydraulic behavior of the substrate. A substrate-specific calibration curve was developed in the laboratory to link the measured soil permittivity to the actual volumetric water content (Kargas et al., 2013;Seyfried et al., 2005). An evaporation pan and a tipping bucket to measure the evaporation demand and the subsurface outflow, respectively, completed the monitoring system connected to an external data logger. The acquisition frequency was set to 1 min.
The assembled testbeds were placed into a climate chamber, where spring-summer Mediterranean climatic conditions were simulated. In particular, air temperatures, relative humidity, and day/night alternations were set to replicate measured data from a weather station in southern Italy (Fig. 2). No radiation intensity control was available in the chamber; therefore, only day/night alternation was simulated. The

Fig. 1.
A schematic of the experimental setup showing the testbeds' structure, GRs' layering, and sensors: (1) a coaxial impedance dielectric reflectometer for the measurement of the volumetric water content in the middle of the soil substrate (z = − 6 cm), (2) thermocouples to measure the temperature in the substrate (z = − 6 cm) and at the bottom of the GRs, (3) a tipping-bucket rain gauge for the measurement of outflow, which is then sampled for tracer and nitrogen analyses, and (4) an evaporation pan to estimate the evaporation demand in the chamber. experiment lasted 41 days and encompassed the monitoring of water contents, subsurface outflow, and the tracer application. Nitrogen compounds in subsurface outflow were analyzed to characterize water flow and nitrogen turnover in the GRs. Two testbeds were irrigated with a 1:2 primary wastewater effluent (1 part pure water + 1 part wastewater), while pure water was applied to the control box. The effluent was characterized by negligible amounts of organic nitrogen and nitrates and an average ammonium concentration of 25 mg/L, which is in line with other studies (e.g., Do et al., 2019;Tao et al., 2014).

Transient water flow and tracer experiments
The irrigation schedule was designed to 1) simulate the real operating conditions of GRs for domestic wastewater treatment, 2) induce an alternation of wet and dry soil conditions, and 3) minimize the plant water stress. In particular, each box was irrigated four times per week for one hour, during which 3 L of applied solution were distributed on the soil surface. The amount of applied water was doubled during the tracer experiments to speed up the tracer's recovery. As shown in Brunetti et al. (2018), such an irrigation pattern can significantly enhance the evaporative cooling efficiency of GRs.
Measured volumetric water contents and subsurface outflow were used to assess the GR's hydrological behavior, particularly its retention capacity and unsaturated conditions in the growing medium. The two tracer experiments were conducted to characterize water flow and solute transport at the beginning and the end of the experimental campaign. Breakthrough concentration curves were reconstructed for each tracer application by measuring solute concentrations in the outflow over time. The first tracer experiment was conducted on June 24th, 2019. A known amount of deuterated water (1 mL, 90atom% deuterated water) and salt (30 g NaCl) was applied to all testbeds together with a mix of wastewater and pure water in a 1:1 proportion (6 L). The same irrigation pattern (6 L per day) was applied in the following days until the outflow salt concentration reached the background concentration. The first tracer experiment lasted 4 irrigation cycles (4 days and 24 L in each box). The second tracer experiment, also lasting 4 irrigation cycles (days), was conducted using the same procedure 5 weeks later (on July 29th, 2019).
The outflow collection in both experiments was done using 50 mL plastic vials at equal intervals (15 min) until the outflow was negligible. The salt concentration was measured using an electrical conductivity (EC) probe. Afterward, an initial salt concentration curve was constructed for each box. Fifty subsamples for each box were taken for the isotopic analysis of deuterated water. Samples were analyzed with laser spectrometry (L 2140-i, Picarro Inc., Santa Clara, CA, USA). To diminish the memory effect and attain more precise measurements, the analysis was set to 11 subsequent applications, while the last three were averaged during the post-processing analysis to get the mean sample concentration. Results given in delta notation were converted into concentrations.

Characterization of the reactive nitrogen transport
Inorganic nitrogen compounds (ammonium, nitrite, nitrate) were monitored weekly in the three boxes' outflow. To ensure the data reliability, sampling was done randomly relative to the irrigation cycles two times per week with 2 to 3 daily samples. The sampling and laboratory analysis were done in replicates to minimize errors and avoid outliers. Plastic vials with 50 mL volume used for sampling were stored in the fridge for a maximum of 3 days due to the unstable nature of nitrogen species.
UV-VIS spectrophotometry was used to measure inorganic nitrogen species. All samples were filtered through 0.1 μm Savana® Lab Disc membrane filters and stored for up to 2 days at 4 • C without acid preservation to reduce the potential interference of the dissolved organic matter. Ammonium and nitrite were measured colorimetrically. The determination of ammonium is based on the reaction with phenol and hypochlorite to give an indophenol blue in an alkaline medium (Berthelot reaction). After the sample preparation, a 30 min reaction time is necessary before measuring the absorbance at 667 nm. The nitrite determination is based on the nitrite's reaction in the acid solution of a primary amine (Sulphanilamide) to form diazonium salt, which is coupled to an aromatic amine to produce the red azo dye, whose absorbance can be measured at 540 nm. The second-derivative ultraviolet spectrophotometric method was used to measure nitrate in the outflow (Ferree and Shannon, 2001). All samples and standards were treated in the same way: 5 mL of sample (or standard) was transferred into a glass vial and mixed thoroughly with 1 mL of sulfuric acid (H 2 SO 4 ) to reduce the pH. Subsequently, each sample's absorbance (or standard) was measured in the range of 200 and 230 nm in the spectrophotometer, and the maximum second derivative was calculated. The calibration curve was fitted to the maximum second derivative of the standard samples and used to calculate outflow nitrate concentrations.

Modeling theory 2.2.1. Water flow and reactive solute transport
Considering the small dimensions of the testbed (Fig. 1), the transport processes are assumed one-dimensional and strictly vertical. Variably-saturated water flow is simulated by solving the Richards equation. Brunetti et al. (2020a) compared multiple hydraulic parameterizations of the growing medium and demonstrated that the unimodal van Genuchten-Mualem (VGM) model (van Genuchten, 1980) the appropriate level of model complexity for the numerical analysis of GRs. Therefore, this parameterization is used in the present study.
Root water uptake is simulated using the approach of Feddes et al. (1978). First, the measured potential evapotranspiration demand in the chamber is partitioned into potential evaporation and transpiration fluxes using the Leaf Area Index (LAI) (Sutanto et al., 2012). LAI is not constant but varies over time due to plant growth, which is assumed to be logistic. The calculated potential transpiration flux is then converted into actual root water uptake by combining the piecewise linear water stress response model proposed by Feddes (Feddes et al., 1978) ) and a root density function that accounts for the root density and growth. A thorough description of the evapotranspiration partitioning and root water uptake is provided in Brunetti et al. (2019).
The reactive transport of nitrogen in the GR is modeled as a sequential reaction chain, in which organic nitrogen (N org ) is first converted into ammonium (NH 4 + ) and then directly oxidized to nitrate (NO 3 -). The intermediate conversion from ammonium to nitrite and from nitrite to nitrate is assumed to be instantaneous and thus neglected in the model. First-(µ) and zero-order (γ) coefficients describe nitrogen turnover and potential production of organic nitrogen by the degradation of plant residues, respectively. The effect of denitrification has been neglected due to the strongly unsaturated conditions in the substrate during the experiment. The transport of the kth solute species is described using the advection-dispersion-reaction equation, assuming that solutes can exist only in the solid and liquid phases. Linear adsorption to the solid phase is considered for N org and NH 4 + (Gärdenäs et al., 2005;Li et al., 2015;Mariano et al., 2016), while 2 H and NO 3 exist only in the liquid phase and cannot react to form other compounds. Unlimited passive uptake of NH 4 + and NO 3 in the liquid phase is considered (Li et al., 2015). The effect of molecular diffusion is neglected.

Numerical domain and boundary conditions
The finite-element model HYDRUS-1D (Šimůnek et al., 2016) is used to solve the Richards and advection-dispersion-reaction equations. HYDRUS has been used successfully to simulate the hydrological performance of multiple Low Impact Development techniques (e.g., Brunetti et al., 2016b;Hilten et al., 2008;Palla et al., 2009). As in Brunetti et al. (2020a), the main modeling assumptions are that the hydraulic effect of the drainage layer is negligible and that water flow and solute transport are strictly vertical. The vertical domain is discretized in 100 one-dimensional elements, which are refined near the soil surface. An "atmospheric" boundary condition is specified on the top boundary, while a "seepage face" is used for the bottom boundary. A spin-up period of 40 days is used to relax the initial conditions. The initial soil nitrogen content is set to 4.1%, according to laboratory measurements. The plant growth is assumed logistic and characterized by a growth rate of 0.15 1/ day (Brunetti et al., 2020a). The initial and final LAIs are set to 1.0 and 3.0, respectively. The simulated root depth varies between an initial measured root depth of 4 cm and a final measured value of 12 cm.

Inverse modeling scenarios
The calibration framework encompasses three modeling scenarios focused on the inverse estimation of soil hydraulic and reactive transport parameters of the three testbeds: control (CR), non-vegetated (NV), and vegetated (VG) (Fig. 1. Measured volumetric water contents, subsurface outflow, tracer breakthrough curve from the first experiment (Fig. 2, and outflow nitrate concentrations are used for model calibration. Outflow ammonium concentrations were negligible; therefore, they were not included in the analysis. The soil hydraulic parameters (θ s , α, n, K s , l) and the solute dispersivity (λ) of the control box were already estimated in Brunetti et al. (2020a) with limited uncertainty. Thus, in the present study and this scenario, these parameters are set to their median value, and the analysis is restricted to the estimation of reactive parameters. Both the root growth and solute and water uptake are neglected in the NV scenario, with potential evapotranspiration converted into potential evaporation and actual evaporation depending on the surface moisture availability. The production of nitrogen from plants' tissue degradation is considered only for scenarios CR and VG.

MULTINEST algorithm
A Bayesian framework is used to calibrate the model and assess its predictive uncertainty. In particular, the nested sampling estimator MULTINEST (Feroz et al., 2009) is combined with experimental observations to inversely estimate the posterior distributions of soil hydraulic and reactive transport parameters for the three scenarios previously described. For a thorough description of the MULTINEST algorithm, please refer to Feroz et al. (2009) andBrunetti et al. (2020b).
The main advantage of nested sampling algorithms is their computational efficiency and the possibility to simultaneously calculate the parameters' posterior distribution and the marginal likelihood, which provide a valuable statistical basis for model selection. Based on a preliminary sensitivity analysis of the algorithm's stability, the number of live points in the MULTINEST was set to 100, and the maximum number of algorithm iterations was set to 200,000. This threshold value was never reached during the analysis, which was performed three times for each modeling scenario to assess the results' accuracy. The stopping criterion is then based on the ratio between the estimated total evidence Z est and the current evidence Z i : (1)

Likelihood function and prior distributions
Error residuals were assumed to be uncorrelated and normally distributed with a constant variance, σ 2 , thus leading to the following log-likelihood function ℓ(Ω): where k is the number of measurements, and H i (Ω) and ỹ i are the ith model realization and its corresponding measured value, respectively. In the Bayesian analysis, the log-likelihood function to be maximized is the aggregated sum of single: where ℓ VWC (Ω) is the log-likelihood of volumetric water content, ℓ OUT (Ω) is the log-likelihood of subsurface outflow, ℓ 2H (Ω) is the loglikelihood of outflow deuterium concentration, and ℓ NO3-(Ω) is the log-likelihood of outflow nitrate concentration. A very low ℓ(Ω) is attributed to both non-convergent model runs and simulations affected by high mass balance errors. The standard deviation in eq. (2) is inferred by the Bayesian framework for each measurement type. Uniform prior distributions are used in the Bayesian analysis. Parameter bounds are provided in Table 1 and are set based on preliminary analyses and a careful literature review (Brunetti et al., 2020a;Brunetti et al., 2020b;Brunetti et al., 2016a;Hanson et al., 2006;Hill et al., 2019;Sandoval et al., 2017). The residual water content, θ r , is set to 0.0, which is a realistic approximation for coarse-textured soils. A high value of the partition coefficient is used for the organic nitrogen to simulate its strong sorption to the solid phase (Brunetti et al., 2020b), which hinders its transport and occurrence in the liquid phase. Thus, only the N org sorbed fraction can react to form NH 4 + . The upper boundary for the saturated water content is set to be lower than the measured porosity to account for the effect of air entrapment (Snehota et al., 2015). Finally, posterior predictive checks are used to "look for systematic discrepancies between real and simulated data" (Gelman et al., 2004). Despite the questionable criticism of using the data twice, we want to emphasize that the posterior predictive checks are a fundamental tool for the analyst to assess the model adequacy (Gelman et al., 2004;Kruschke, 2013).

Nitrate leaching from Green Roofs: Global sensitivity analysis
After the model calibration and the uncertainty assessment, a Global Sensitivity Analysis is performed to identify the most important factors driving the leaching of nitrates from the GR, which is directly linked to the system's environmental impact. To include the vegetation's effect, scenario VG is used as the basis for the numerical analysis. This choice is intended to provide a statistical basis to enlighten the relative importance of different physicochemical processes on the production and transport of nitrates in the GR.
The Random Balance Designs Fourier Amplitude Sensitivity Test (RBD-FAST) Tissot and Prieur, 2012) is coupled with HYDRUS. The method combines the accuracy of the classic Fourier Amplitude Sensitivity Test (Saltelli et al., 1999) with the computational efficiency of Satterthwaite's random balance designs (Satterthwaite, 1959) to estimate the first-order effect (S1) in a variance-based context (Saltelli et al., 2006). Higher S1 values are attributed to more influential parameters. The main advantage of the RBD-FAST method is that the total number of model runs is reduced down to N instead of d*N (like in Sobol' or FAST methods), where d is the number of factors investigated.
In the present study, a convergence analysis is used to determine N. Nonconvergent runs are excluded a posteriori from the computation of the sensitivity indices. Parameter bounds are the same as used for the Bayesian analysis and reported in Table 1.

Water flow and solute transport
Results of the model calibration are shown in Fig. 3, which compares the observations and the posterior predictive checks (grey lines) of simulated subsurface outflow [(A) and (D) in Fig. 3 Fig. 3 Fig. 3] for all scenarios. Model predictions are obtained by random sampling of 100 solutions from the posterior parameters' distributions, whose quantiles are listed in Table 1.

], and volumetric water contents [(C) and (F) in Fig. 3] for scenarios NV and VG, and outflow nitrate concentrations [(G), (H), and (I) in
The fitting quality is satisfactory, and the predictive uncertainty is limited for both scenarios, although the model overestimates measured subsurface outflow from the non-vegetated testbed. This deviation accumulates at the beginning of the numerical simulation. Still, it doesn't grow significantly over time, suggesting a potential overestimation of the inflow during the tracer application at the beginning of the experimental campaign. Conversely, the quality of the fitting of the volumetric water content slightly degrades over time, with the model underestimating the observations. Inaccuracies in simulated actual evaporation can explain this deviation, whose magnitude remains anyway limited.
Simulated actual evapotranspiration is also partially responsible for the minor differences between simulated and observed volumetric water Table 1 Model parameters, their bounds, and 5%, 50%, and 95% quantiles of the parameters' posterior distributions calculated from the Bayesian analysis. The last column reports the first-order sensitivity indices, S1, calculated using the RBD-FAST method to assess the influence of different factors on the leaching of nitrates from the vegetated GR.  - † The soil hydraulic parameters and solute dispersivity for this scenario are set according to Brunetti et al. (2020a).
contents in scenario VG, mainly encountered after irrigation events. Nevertheless, the fitting quality remains satisfactory, suggesting that the calibrated model reproduces the hydraulic behavior of the vegetated GR well. The accurate description of subsurface outflow further confirms this. In contrast, the tracer breakthrough curve is slightly overestimated, as reflected in the wider confidence interval for the solute dispersivity (Table 1. This deviaiton can also result from mass differences when injecting the tracer cocktail or by inhomogeneous distribution of the irrigation water on the soil surface. The estimated soil hydraulic parameters are characterized by narrow confidence intervals for both inverse modeling scenarios, NV and VG, which exhibit similar estimated saturated water contents, θ s . These values are both considerably lower than that reported by Brunetti et al. (2020a) for the control box. Significant differences also emerge for the saturated hydraulic conductivity and the tortuosity. Altogether, these deviations can be partially explained by soil heterogeneity, air entrapment, and nonuniform application of irrigation water on the soil surface.
The estimated shape parameters, α and n, indicate a soil characterized by a sharp air-entry point and a moderate retention capacity, typical of GR substrates (Brunetti et al., 2016a;Hill et al., 2019;Sandoval et al., 2017). The calibrated Feddes' parameters reflect the nearly optimal transpiration rates of the GR vegetation under wet soil conditions and agree with previous studies (Brunetti et al., 2016a;Kuronuma and Watanabe, 2017).

Reactive nitrogen transport
Nitrogen turnover is well predicted in all inverse modeling scenarios. The highest predictive uncertainty is encountered for scenario NV [(H) in Fig. 3], while the scenarios CR [(G) in Fig. 3] and VG [(I) in Fig. 3] exhibit both narrow prediction bands. This trend was already observed for water flow and solute transport, suggesting experimental inaccuracies in the non-vegetated testbed. Nevertheless, the description of the nitrogen reactive transport remains sufficiently accurate also for this scenario.
Comparing the estimated parameters posterior distributions holds more interest and clarifies the role of different physical and biological processes. First, the very low values estimated for γ Norg indicate that the release of nitrogen from the decomposition of plant residues was limited during the experiment. Although this was expected due to the limited duration of the experimental campaign, this process cannot be excluded in the long term and under field operating conditions, where the seasonality can trigger changes in the plant status. Interestingly, the estimated degradation rate of organic nitrogen, µ SNorg , for the vegetated box is significantly higher than that calculated for the control box. This suggests that wastewater application induces optimal conditions for the microbial conversion of N org to NH 4 + (Ibekwe et al., 2018). The NH 4 + Fig. 3. Model calibration results. (A) and (D): measured (red lines) and modeled (grey lines) cumulative subsurface outflow for the NV and VG scenarios, respectively; (B) and (E): measured (blue circles) and modeled (grey lines) deuterium concentrations in the outflow for the NV and VG scenarios, respectively; (C) and (F): measured (black lines) and modeled (grey lines) volumetric water contents at z = -6 cm for the NV and VG scenarios, respectively; (G), (H), and (I): measured (black circles) and modeled/grey lines) nitrate concentrations in the outflow for the CR, NV, and VG scenarios, respectively. Model predictions are obtained by random sampling of 100 solutions from the posterior parameters' distributions. (G), (H), and (I) have the same scale to visually appreciate the difference in nitrates leaching between testbeds. RMSE is the root mean square error calculated from the median solution of the posterior parameters distribution.
partition coefficient cannot be inferred for the CR and NV scenarios. Still, its estimated value for the VG scenario indicates strong sorption of this compound to the solid phase, compatible with the high organic matter content of the substrate (Wang and Alva, 2000). The high to moderate uncertainty of this estimated parameter is mainly a consequence of the lack of information about the ammonium distribution in the soil profile. The NH 4 + degradation rate in the liquid phase, µ LNH4+ , exhibits high uncertainty and cannot be inferred from the observations. Conversely, the high value estimated for µ SNH4+ , which, together with µ SNorg , indicates high nitrification processes in the vegetated testbed. The comparison between the estimated standard deviations of nitrate measurements, σ NO3-, confirms experimental inaccuracies for the nonvegetated testbed and might cause the high parameter uncertainty for this scenario. Table 2 reports the cumulative N fluxes in the GRs for different scenarios calculated by averaging the results of 100 numerical simulations obtained by random sampling from the posterior parameters' distributions. The results indicate that the N leaching potential is significantly higher in the non-vegetated testbed, where it reaches 94% of the total injected N, compared to the 67% estimated for the vegetated box. In the latter scenario, the roots take up 32% of the injected N, while N input from decomposition plays a minor role. Conversely, this process is more important for the control scenario, where 63 and 30% of the generated N are leached and taken up by roots, respectively. These percentages are similar to those obtained for the vegetated scenario and suggest that root solute uptake can mitigate N leaching from GRs.

Key processes in nitrogen turnover
The comparison between NO 3 outflow concentrations and the model calibration results raises interesting discussion points on the role of different processes on the nitrogen reactive transport: • The control box (irrigated with tap water) is a source of nitrates. The average NO 3 outflow concentration is 3.6 ± 0.94 mg/L, which is above the reported threshold concentration of 2.6 mg/L that accelerates eutrophication in water bodies (Bartley and Johnson, 2006). This behavior was already observed by Vijayaraghavan et al. (2012), who carried out an extensive experimental campaign on the assessment of GRs water quality. The authors reported similar concentration values when GRs were irrigated solely with tap water and concluded that the nitrate release is mainly due to the degradation of the soil substrate's organic compounds. This agrees with the results of the present study that shows a decreasing trend in nitrate concentrations [(G) in Fig. 3] caused by the consumption of the initial soil organic nitrogen, which is not replenished by the degradation of plant residuals (low γ Norg in Table 1. This behavior was also observed by Monterusso et al. (2004) and suggested the need to periodically fertilize traditional GRs to restore nutrient levels and support plants' growth. Using a numerical model can help optimize fertilization to avoid excess nutrients leaching (Berndtsson et al., 2006). Finally, the narrow estimated uncertainty for µ SNorg confirms that organic nitrogen turnover is the dominant process in the release of nitrates from the control box. • The non-vegetated GR releases the highest amount of nitrates, suggesting an important role of root solute uptake for the nitrogen removal from GRs irrigated with domestic wastewater. This uptake is favored by moderately wet conditions of the soil substrate during the experiment [(F) in Fig. 3]. Vijayaraghavan et al. (2012) reported similar differences between vegetated and non-vegetated extensive GRs. Similarly, Berndtsson et al. (2006) highlighted the role of vegetation on nitrogen retention in extensive GRs. NO 3 leaching from the non-vegetated testbed exhibits a decreasing trend [(G) in Fig. 3], which is again compatible with gradual degradation of the initially present soil organic nitrogen. This trend is not observed for the vegetated box and can be partially explained by plant uptake, which hampers the initial accumulation of nitrates. Afterward, the applied wastewater replenishes the soil nitrogen and induces optimal nitrification conditions (high µ SNH4+ in Table 1, thus maintaining an approximately constant outflow concentration.

Global Sensitivity analysis
The previous discussion on the experimental data and model calibration outcomes identified root solute uptake and soil organic compounds' degradation as important processes in nitrogen turnover in the GRs. To further confirm or disprove these findings, a Global Sensitivity Analysis is performed to determine the factors driving nitrate leaching from the vegetated GRs. The convergence analysis indicated that 3000 model executions were sufficient to obtain stable estimations of the firstorder sensitivity indices (reported in the last column of Table 1 for all investigated parameters. Seven factors out of 15 exhibit an appreciable influence (S1 >0.1%) on nitrate leaching, but only five are characterized by S1 >1%: γ Norg , POpt, µ SNH4+ , µ SNorg , and P0 (in order of importance). The soil hydraulic and solute non-reactive transport parameters have negligible influence in the investigated range (Table 1. This may be partially due to the limited depth (i.e., 12 cm) of the soil substrate, which masks the effect of transport processes, and the controlled laboratory conditions, which reduce heterogeneity. Similar reasoning can be applied to the ammonium distribution coefficient, whose retardation effect on the nitrate transport cannot be fully appreciated in a shallow soil profile.
As expected, the first-order degradation coefficients, µ SNH4+ and µ SNorg , have an appreciable impact on nitrate leaching. However, their effect is topped by γ Norg and POpt, which are the two most influential factors. These results confirm that root solute uptake and soil organic matter degradation are important processes in the GR nitrogen balance. To further clarify their functioning, a local sensitivity analysis is performed using the median solution in Table 1 as a reference and changing alternatively only γ Norg and POpt. Fig. 4 shows the resulting cumulative root  Fig. 4]. The comparison confirms the bigger influence of γ Norg . High production of organic nitrogen is positively correlated with an increase in the soil's nitrate level, which cannot be completely offset by an increase in root solute uptake. The latter depends on the plant uptake capacity, which is regulated by the transpiration demand, the vegetation's characteristics, and the soil volumetric water content in the root zone (Brunetti et al., 2021). For instance, very wet conditions can induce waterlogging effects and drastically reduce root water and solute uptake. The Feddes' parameter POpt is a threshold pressure head value, above which actual transpiration drops due to anaerobic conditions. The analysis shows that when POpt is reduced [(C) in Fig. 4], nitrate uptake declines due to the wet conditions observed during the experiment [(F) in Fig. 3]. As a consequence, nitrate leaching is higher than in other modeling scenarios.
To summarize, the main findings of the Global Sensitivity Analysis are: • Organic nitrogen production due to exogenous variables further increases nutrient leaching from the GR, reducing its environmental benefits. Therefore, installing plants characterized by a minimal potential of tissue decomposition and high resistance and designing an irrigation system to avoid extended periods of water stress is preferable. At the same time, fertilization or soil amendments should be carefully evaluated and potentially excluded when domestic wastewater is applied. When this is not possible for various reasons, the fertilizer application should be properly designed to minimize nutrients leaching (Clark and Zheng, 2014) by taking into account the in-situ conditions of the GR. The present study demonstrates that numerical models are a valuable and recommended tool for this task. • Root solute uptake enhances the nitrogen removal efficiency from domestic wastewater. Thus, the system should be designed to maximize transpiration by reducing water stress (both due to waterlogging or dry conditions) and supporting plants' health. Again, we advocate using a numerical model, which can provide a comprehensive overview of physical processes in the GR.

Strengths and limitations
We consider this study to be innovative in mainly two aspects: • The numerical analysis of nitrogen turnover in GRs irrigated with domestic wastewater is a novel application in this research field. The joint use of Bayesian inference techniques and sensitivity analyses provides a statistical basis for uncertainty assessment and ranking of the effects of multiple physical processes on nutrients leaching from the GRs, respectively. • The model calibration is performed on a comprehensive dataset from a controlled laboratory experiment on multiple testbeds, which included not only nitrogen observations but also measurements of volumetric water contents, subsurface outflow, tracer experiments, and atmospheric variables in a climate chamber. This allows to narrow down uncertainties and better understand the physical functioning of the system.
Still, the study also possesses some limitations and open questions to be solved in the future. First, the experimental campaign's duration was too limited to observe long-term processes such as the alteration of the soil hydraulic properties due to biofilm formation or wastewater toxicity effects on the vegetation, which would change both the hydraulic regime and the root uptake capacity. Second, it was impossible to precisely infer certain model parameters (e.g., ammonium partition coefficient) due to both operating conditions and experimental inaccuracies and limitations. Additional measurements of ammonium concentrations in the soil profile could potentially mitigate this problem. Finally, the GRs' morphology, simulated atmospheric conditions, and irrigation pattern are not representative of a broad class of field-scale operating conditions, and more tests under variable atmospheric conditions can give more insights into the physical, chemical, and biological dynamics of GRs, such as denitrification processes which might occur during the long-duration and low-intensity rainfall events (Qiu et al., 2020) or due to water-soluble organic carbon from plant residues during initial decomposition (Surey et al., 2020). Furthermore, precipitation can affect the nitrogen turnover by increasing the nitrate leaching potential and soil water content, which influences the microbial degradation processes (e.g., Maenhout et al., 2018). However, we want to emphasize that observations' information content was sufficient to estimate multiple parameters with limited uncertainty (Table 1.

Conclusions
This study aimed to expand the current knowledge on the potential use of Green Roofs for domestic wastewater treatment, with particular emphasis on the nitrogen reactive transport and nitrate leaching. A Bayesian framework was used to calibrate a mechanistic numerical model against multiple observations from a controlled laboratory experiment. This analysis was complemented by a Global Sensitivity  Table 1 was used as a reference, and only the values of γ Norg and K f were increased 25 and 50 times to have a fair comparison between the two factors.
Analysis to further enlighten the role of multiple physical processes. Results indicate that GRs quickly nitrify wastewater ammonium, and resulting nitrate leaching is negatively correlated with the vegetation's root solute uptake capacity. Maximization of this effect is desirable and could be achieved by selecting appropriate plant species and reducing water stress periods with a proper irrigation schedule. This would simultaneously support vegetation's health and limit the decay and decomposition of plants' tissues, which the analysis showed can potentially influence the nitrogen balance in the GRs. Numerical models are certainly a valuable tool for these tasks, and to understand and design GRs for wastewater treatment. Based on these results, we recommend extending the experimental campaign in future studies to observe the long-term behavior of the system under different operating conditions. We further recommend investigating the potential use of GRs for the decentralized removal of other compounds, such as pharmaceuticals, which show significant persistence and occurrence in the environment (Brunetti et al., 2022). GR vegetations could potentially uptake and metabolize many of these neutral compounds that can easily cross the root membrane (Brunetti et al., 2021).

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.