Chemical Engineering Journal 435 (2022) 134733 Contents lists available at ScienceDirect Chemical Engineering Journal journal homepage: www.elsevier.com/locate/cej Systematic calibration of N2O emissions from a full-scale WWTP including a tracer test and a global sensitivity approach Borja Solís a, Albert Guisasola a, *, Maite Pijuan b, c, Lluís Corominas b, c, Juan Antonio Baeza a a GENOCOV. Departament d’Enginyeria Química, Biològica i Ambiental. Escola d’Enginyeria. Universitat Autònoma de Barcelona, 08193, Bellaterra (Barcelona), Spain b Catalan Institute for Water Research (ICRA), Scientific and Technological Park of the University of Girona, Emili Grahit 101, 17003, Girona, Spain c Universitat de Girona, Plaça de Sant Domènec 3, 17004 Girona, Spain A R T I C L E I N F O A B S T R A C T Keywords: Nitrous oxide (N2O) is a greenhouse gas (GHG) emitted during biological nitrogen removal from wastewater Nitrous oxide treatment plants (WWTPs). Some modelling tools have been proposed to predict N2O emissions during the design WWTP modelling and operation of WWTPs. In this study, the novel ASM2d-N2O model, which accounts for the production of N2O Global sensitivity analysis in nutrient removal WWTPs, was used to study the associated emissions from a full-scale WWTP with two in­ Dynamic calibration dependent lines. Firstly, the hydraulics of the WWTP was characterized by a residence time distribution test, showing the flow was equally divided into the two treatment lines (49.3 vs. 50.7%), that each reactor worked as an ideal continuous stirred tank reactor and the secondary settler model flux was similar to a plug-flow reactor. The ASM2d-N2O model was then calibrated using experimental data obtained under dynamic conditions. A global sensitivity analysis was used to select, among 59 model parameters, five candidates that resulted to be related to nitrifying organisms. Different parameter subsets up to four parameters were evaluated, being the subset [µNOB, qAOB_AMO, KO2_NOB, KNO2_NOB] the best, achieving 53.3% reduction of the calibration cost function. The model fit obtained provided a reasonably description of nutrients and N2O emission trends, considering the inherent operational variability suffered in full-scale WWTPs. Finally, a simulation-based study showed that, for the given WWTP and operational conditions, an unbalanced distribution of flow-rate between the two treatment lines did not result in a significant increase on N2O emissions. The results obtained show that this model can be a suitable tool for predicting N2O emissions in full-scale WWTPs, and can therefore be used to find operational conditions that help to minimise these emissions. Abbreviations: AER, Aerobic reactor; ANA, Anaerobic reactor; ANX, Anoxic reactor; AOB, Ammonia Oxidizing Bacteria; AS, Activated Sludge; ASM, Activated Sludge Model; ASM2, Activated Sludge Model number 2; ASMN, Activated Sludge Model for Nitrogen; BNR, Biological Nitrogen Removal; Br− , Bromide; C, Carbon; CCF, Calibration Cost Function; CCFN2O, Calibration Cost Function for Nitrous oxide; CCFNH4, Calibration Cost Function for Ammonium; CCFNO2, Calibration Cost Function for Nitrite; CCFNO3, Calibration Cost Function for Nitrate; CritD, FIM D criteria; CritModE, FIM modified E criteria; COD, Chemical Oxygen Demand; CODfil, Filtered COD; CODfil, eff, Effluent filtered COD; CODsol, Soluble COD; CODtot, Total COD; CSTR, Continuous Stirred Tank Reactor; DO, Dissolved Oxygen; FIM, Fisher Information Matrix; GHG, Greenhouse gas; GSA, Global Sensitivity Analysis; GWP, Global Warming Potential; HCCF, Hydraulic Calibration Cost Function; HD pathway, Heterotrophic Denitrification pathway; IR, Internal Recirculation; IWA, International Water Association; KBr, Potassium Bromide; MAE, Mean Absolute Error; MC, Monte Carlo; N, Nitrogen; N2, Nitrogen gas; N2O, Nitrous oxide; N2O-EF, Nitrous oxide emission factor; N2O-N, Nitrous oxide nitrogen; ND pathway, Nitrifier Denitrification pathway; NH2OH, Hydroxylamine; NH4+, Ammonium; NN pathway, Nitrifier Nitrification pathway; NO, Nitric oxide; NO2–, Nitrite; NO3–, Nitrate; NOB, Nitrite Oxidizing Bacteria; NormD, FIM normalized D criteria; P, Phosphorus; PAO, Polyphosphate Accumulating Organisms; PCCF, Preliminary Calibration Cost Function; PO43− , Phosphate; R2, Coefficient of determination; RDE, FIM ratio between NormD and CritModE; RMSE, Root Mean Square Error; RSA, Regional Sensitivity Analysis; RTD, Residence Time Distribution; SI, Supplementary Information; Si, Sensitivity index; SS, Steady-State; TKN, Total Kjeldahl Nitrogen; WWTP, Wastewater Treatment Plant; XMeOH, Ferric Hydroxide state variable. * Corresponding author. E-mail addresses: Borja.solis@uab.cat (B. Solís), albert.guisasola@uab.cat (A. Guisasola), mpijuan@icra.cat (M. Pijuan), lcorominas@icra.cat (L. Corominas), JuanAntonio.Baeza@uab.cat (J.A. Baeza). https://doi.org/10.1016/j.cej.2022.134733 Received 28 May 2021; Received in revised form 26 December 2021; Accepted 12 January 2022 Available online 18 January 2022 1385-8947/© 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). B. Solís et al. Chemical Engineering Journal 435 (2022) 134733 1. Introduction Hydraulics is considered relevant from two perspectives that affect the performance of the plant: pollutant degradation depends on the flow Over the past years, concerns regarding the sustainability of waste­ rate and the control of the process is based on concentrations measured water treatment plants (WWTPs) have increased, with particular in specific points in the reactor [17]. Classical tools as residence time emphasis on the carbon (C) footprint. Nitrous oxide (N2O) is a green­ distribution (RTD) techniques with tracer experiments have been house gas (GHG) that is emitted during biological nitrogen removal applied to investigate the behaviour of full-scale WWTPs, in order to (BNR) in WWTPs. Due to its high global warming potential, 265 times detect hydraulic problems such as plug-flow behaviour [18], unopti­ higher than that of carbon dioxide (CO2) [1], the C footprint of WWTPs mized inlet flow distribution ratio [19], short-circuiting flows and dead are highly sensitive to N2O emissions [2]. In addition, N2O is an ozone zones [20,21]. Studies in full-scale WWTPs have demonstrated that layer depletion gas [3]. Measurement campaigns on full scale WWTPs these hydrodynamic anomalies can reduce pollutant removal and in­ have shown high variability on the measured N2O emissions, with a N2O crease management costs [21]. Moreover, biased sampling strategies not emission factor (N2O-EF, defined as the fraction of influent nitrogen load considering hydraulic RTD and daily flow variations and concentrations emitted as N2O) ranging between 0.01% and 1.8%, and in some cases can also lead to erroneous conclusions [22]. In the present study, where even higher than 10% [4–6]. the plant consists of two parallel lines, there was concern about the The biological pathways for N2O production during BNR are related possible effect on N2O emissions if there was an imbalance between the to nitrification and denitrification processes [4], being the first one the flows treated in both lines. Plant overloading or changes in plant load major contributor to N2O emissions in full–scale plants. Nitrification have been one of the operating conditions believed to have a significant consists of two coupled processes: i) ammonium (NH4+) oxidation to impact on N2O emissions [16], and hence, a maldistribution of flow hydroxylamine (NH2OH), nitric oxide (NO) and nitrite (NO2–) by rates between the lines could cause an overloaded line significantly ammonia oxidizing bacteria (AOB) [7] through the nitritation process increasing its N2O emissions. and ii) nitrite (NO2–) oxidation to nitrate (NO3–) by nitrite oxidizing Therefore, the main objective of this study was to comprehensively bacteria (NOB) through the nitratation process. AOB are responsible of calibrate the ASM2d-N2O model using dynamic data from a full-scale N2O production during BNR by two possible pathways. On the one hand, municipal WWTP. Firstly, an RTD experiment was performed to cali­ N2O can be produced through the biological NO reduction due to brate WWTP hydraulics, as an essential initial step to obtain accurate incomplete NH2OH oxidation: the NH2OH incomplete oxidation predictions of the N-species and N2O emissions. Then, an experimental pathway (NN pathway) [4]. On the other hand, AOB can produce N2O campaign was carried out during three days in different zones of the by the sequential reduction of NO2– to NO and N2O: the nitrifier deni­ WWTP to capture its N2O emission dynamics. The parameter subset trification pathway (ND pathway). There exists a third biological N2O selection to calibrate the ASM2d-N2O model was obtained through a production pathway during the denitrification process (HD pathway). global sensitivity analysis (GSA) method, to capture the interactions N2O is formed as an obligate intermediate of the four step reduction between the parameters and to select the parameter subset indepen­ reactions of NO3– to NO2– and then to NO, N2O and finally dinitrogen gas dently from the initial parameter value. The identifiability of the (N2) [4,8,9]. parameter subsets was evaluated using the Fisher information matrix Mathematical models have been widely applied to the prediction of (FIM) methodology. Finally, a simulation-based study was carried out to nitrogen (N) and phosphorus (P) removal in WWTPs since the devel­ assess the possible effect that an incorrect flow distribution between opment of the activated sludge models (ASM) by the International Water treatment lines had on the N2O emissions from the WWTP due to Association (IWA) [10]. Over the past years, modelling BNR has gained possible treatment line overloading and changes in treatment plant more attention in view of a better understanding of N2O production, loads. accumulation and emission. The ability to predict N2O emissions serves as a method for verifying hypothesis related to fundamental mechanisms 2. Materials and methods for N2O production, and it can be used to anticipate N2O emissions in the design and operation of WWTPs, as well as in the design of potential 2.1. Girona WWTP description mitigation strategies [6,11,12]. Different models have been developed aiming at predicting lab-scale The full-scale urban WWTP is located in Girona (Catalonia, Spain) or full-scale N2O emissions [6,12–16]. These models are based on with a design capacity of 206250 population equivalent and 55000 m3 different assumptions, incorporating one, two or three of the biological d− 1. A schematic overview of the WWTP process layout is shown in pathways for N2O production (ND, NN or HD). Regarding nitrification, it Fig. 1. After the pre-treatment (grit and gross removal), the influent is was seen that the models only including a single N2O production distributed in three rectangular lamella-plate primary clarifiers. The pathway could not explain all the experimental data in the literature. biological section has a 5-stage Bardenpho configuration. It consists of Therefore, it was formulated that both NN and ND can occur at the same two main treatment lines with seven separated reactors in each line. The time depending on the operating conditions [14]. Among the different wastewater flows through an anaerobic reactor (ANA1, 1335 m3), an published N2O models, the ASM2d-N2O model developed by Massara anoxic reactor (ANX1, 4554 m3), three aerobic reactors with the same et al. [16] is an ASM type model that: includes N, P and organic matter volume (AER1 to AER3, 1929 m3), a second anoxic reactor (ANX2, 1276 removal; integrates all the microbial pathways for N2O production and m3) and finally a fourth aerobic reactor (AER4, 1409 m3). The internal consumption; contains N2O stripping modelling and estimates the N2O- recycle (IR) flows from the third aerobic reactor (AER3) to the first EF under a wide range of operating conditions. Therefore, the ASM2d- anoxic reactor (ANX1) and it is independent in each line (Fig. 1). In N2O model is a promising tool for developing N2O mitigation strategies addition, sodium aluminate is injected into each IR stream to favour during full-scale WWTP. However, although the different sub-models of chemical P precipitation. The flow is mixed again at the outlet of the last the whole ASM2d-N2O were calibrated separately, the prediction aerobic zone and divided into three parallel secondary clarifiers (5332 capability of the ASM2d-N2O model has not been proved during full- m3 each settler) where the biomass is separated from the treated scale treatment. effluent. The concentrated outflows from each secondary clarifier are With regard to the modelling of full-scale WWTPs, the variability of mixed and, after extracting the purge flow, recycled to the influent of the the operational parameters or the hydraulic conditions of the reactor biological reactors. Aeration is supplied by blowers in the aerobic zones. itself is often not considered, and simplifications of ideality end up being Each line has a blower whose air is sparged through diffusers. used, modelling the reactors as continuous stirred tank reactor (CSTR) with perfectly known inlet and outlet flow rates. However, hydraulics should be considered to obtain more realistic model predictions. 2 B. Solís et al. Chemical Engineering Journal 435 (2022) 134733 Fig. 1. Plant layout of the Girona WWTP (PC = Primary Clarifier, SC = Secondary Clarifier, ANA = Anaerobic reactor, ANX = Anoxic reactor and AER = Aero­ bic reactor. 2.2. Hydraulic characterization procedure 2.3. Experimental data campaign Tracer experiments to characterize the hydraulics of the secondary A three-day data sampling campaign was carried out from 18 to 20 treatment of the WWTP (i.e. biological reactors and secondary clarifiers) July 2017 to calibrate the ASM2d-N2O model at the Girona WWTP were carried out with potassium bromide (KBr). The objectives of the under dynamic conditions. Figure S1 in the supplementary information hydraulic characterization were: 1) to determine the flow distribution (SI) shows a scheme of the WWTP with the locations and data collected between both treatment lines in the biological reactor, as the plant op­ during the experimental data campaign. Only one biological treatment erators suspected that it was not equally distributed in each treatment line was sampled (Line 1, Fig. 1). The data collected are summarized line and 2) to understand the hydraulics of each reactor to identify below: possible dead-zones. An amount of 24.9 kg of KBr (16.7 kg of Br− ) was added to the influent of the primary clarifiers. During the experiment, - Chemical analysis: Grab samples were collected to analyse NH4+, samples were taken from different zones of the two treatment lines of the NO2–, NO3– and PO43− by ion chromatography (ICS5000, DIONEX) biological reactor and in the effluent. The locations of the tracer pulse at different locations and intervals. The rationale of the sampling and the different tracer samples zones are indicated in Figure S1 in frequency was to be able to relate N2O emissions to N-species. section S1 of the Supplementary Information (SI). Therefore, the sampling frequency was increased in the reactors were The model structure used to determine the RTD of the activated N2O emissions were monitored. Grab samples were taken 5 times a sludge (AS) system of the Girona WWTP was the n-tanks in series [23], a day at 3 h intervals (from 9:00 to 18:00) plus midnight for 3 days at widely used method in the calibration of full-scale WWTPs [18,24,25]. reactors ANX1, AER1, AER3 and AER4. At reactors ANA1 and ANX2, The pulse was added in the influent of the primary clarifiers but, for the grab samples were taken less frequently, only 3 times a day. In simulation purposes, the data at the effluent of the primary clarifiers addition, two refrigerated automatic samplers also took samples were used as the input simulation because modelling of the complex every 3 h during the 3 days of the data campaign (8 times per day) at hydraulic lamella-plates clarifiers of the primary settler was out of the the influent of the biological reactors and in AER2 reactor. Br− was scope of this work and was not needed for the WWTP calibration. As an measured using a Dionex DX-500 ion chromatograph (Dionex Corp., initial approach, each vessel was considered to be an ideal CSTR and the Sunnyvale, CA, USA). The method detection limit (MDL) was 0.045 influent flow was equally divided into each treatment line. Hence, each µg/l for Br− . The practical limit of quantification (LOQ = 5 MDL) was line of the secondary treatment was initially simulated as 8 tanks in 0.23 mg/l for Br− . series (7 reactors plus the secondary settler also simulated as a single - Online sensors: NH4+, COD, pH and temperature at the bioreactor tank). The mass balance of the tracer concentration in each tank is inlet were continuously monitored utilizing two on-line ion-selective described by Equation (1): electrodes (ammo::lyserTM) coupled to a monitoring station (S::CAN ( ) Messtechnik GmbH, Austria). In addition, three Dissolved Oxygen dCi Vi ⋅ = Qi ⋅ Cin,i − Ci (1) (DO) probes (oxi::lyserTM) were installed in AER1, AER2 and AER3 dt coupled to the same S::CAN monitoring station. Where: Vi is the volume of reactor “i”; Qi is the volumetric flowrate in - Gas emission measurements: N2O emissions were measured in reactors reactor “i”, Cin,i is the inlet tracer concentration and Ci the tracer con­ AER1, AER2 and AER3 using a system with three gas collection centration in reactor “i”. hoods. The gas collected in the hoods was coupled with a monitoring The calibration cost function of the hydraulic experiment to be unit to log gas temperature, pressure and flowrate. Part of the gas minimized was named HCCF (Equation (2)): was pumped to a conditioning unit (M&C Tech group) to remove humidity and particles and finally the gas was conducted to an online √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ∑ analyser (Horiba VA3000). N2O concentration (in ppmv), pressure, N max|CANA1− L1 | ( )2 HCCF = ⋅ Ci,exp − Ci,model (2) i=1 max|Ci | gas flowrate and temperature were logged at 15 s intervals. The analyser measured only one hood at a time in 20-minute intervals Where: “i” refers to each reactor; Ci,exp is the experimental concentra­ between each reactor. A detailed methodology for the gas emission tion; Ci,model is the model prediction and max|CANA1− L1 |/max|Ci | is a measurement system can be found in Ribera-Guardia et al. [26]. weight factor to normalize the concentrations of each reactor to those of - Data from the WWTP SCADA system: data of the hydraulics (biological ANA1 from the biological treatment line one. Thus, each of the experi­ influent, internal and external recirculation and wastage flowrates) mental inputs has the same influence on the HCCF. and DO of the aerobic reactors were collected from the SCADA sys­ tem of the WWTP. 3 B. Solís et al. Chemical Engineering Journal 435 (2022) 134733 Table 1 Symbols, description, default values at 20 ◦ C, units and variation range of the ASM2d–N2O parameters included in the GSA. Parameter Description Default value at 20 ◦ C Units Min/max range − 1 KH Hydrolysis rate constant 3 d 1.5 / 4.5 KO2_H Saturation/inhibition coefficient for O2 0.2 g O2 m− 3 0.1 / 0.3 Kx_H Saturation coefficient for particulate COD 0.1 g XS (g XH) − 1 0.05 / 0.15 nNO3_H Anoxic hydrolysis reduction factor 0.6 – 0.3 / 0.9 nNO2_H Anoxic hydrolysis reduction factor 0.6 – 0.3 / 0.9 KNO3_H Saturation/inhibition coefficient for NO3– 0.5 g N m− 3 0.25 / 0.75 KNO2_H Saturation/inhibition coefficient for NO2– 0.5 g N m− 3 0.25 / 0.75 nfe_H Anaerobic hydrolysis reduction factor 0.4 – 0.2 / 0.6 µH Maximum growth rate on substrate 6 g XS (g XH) − 1d− 1 3/9 KO2 Saturation/inhibition coefficient for O2 0.1 g O2 m− 3 0.05 / 0.15 KF Saturation coefficient for growth on SF 20 g COD m− 3 10 / 30 KNH4 Saturation coefficient for NН4+ (nutrient) 0.05 g N m− 3 0.025 / 0.075 KP Saturation coefficient for PO43− (nutrient) 0.01 g P m− 3 0.005 / 0.015 KALK Saturation coefficient for alkalinity (HCO3–) 0.1 mole HCO3– m− 3 0.05 / 0.15 KA Saturation coefficient for growth on acetate SA 20 g COD m− 3 10 / 30 KNO3 Saturation/inhibition coefficient for NO3– 0.5 g N m− 3 0.25 / 0.75 KNO2 Saturation/inhibition coefficient for NO2– 0.5 g N m− 3 0.25 / 0.75 nNO3_D Reduction factor for denitrification 0.28 – 0.14 / 0.42 qfe Maximum rate for fermentation 3 g SF (g XH) − 1d− 1 1.5 / 4.5 Kfe_H Saturation coefficient for fermentation of SF 4 g COD m− 3 2/6 bH Rate constant for lysis and decay 0.4 d− 1 0.2 / 0.6 nG3 Anoxic growth factor (NO2–→NO) 0.16 – 0.08 / 0.24 nG4 Anoxic growth factor (NO → N2O) 0.35 – 0.175 / 0.525 nG5 Anoxic growth factor (N2O → N2) 0.35 – 0.175 / 0.525 KS3 Half-saturation coefficient for substrate 20 g COD m− 3 10 / 30 KS4 Half-saturation coefficient for substrate 20 g COD m− 3 10 / 30 KS5 Half-saturation coefficient for substrate 40 g COD m− 3 20 / 60 KNO2_Den Half-saturation coefficient for NO2– 0.2 g N m− 3 0.1 / 0.3 KOH4 Half-saturation coefficient for O2 0.1 g O2 m− 3 0.05 / 0.15 KN2O_Den Half-saturation coefficient for N2O 0.05 g N m− 3 0.025 / 0.075 KOH3 Half-saturation coefficient for O2 0.1 g O2 m− 3 0.05 / 0.15 KNO_Den Half-saturation coefficient for NO 0.05 g N m− 3 0.025 / 0.075 KOH5 Half-saturation coefficient for O2 0.1 g O2 m− 3 0.05 / 0.15 KI3NO NO inhibition coefficient (NO2–→NO) 0.5 g N m− 3 0.25 / 0.75 KI4NO NO inhibition coefficient (NO → N2O) 0.3 g N m− 3 0.15 / 0.45 KI5NO NO inhibition coefficient (N2O → N2) 0.075 g N m− 3 0.038 / 0.112 µAOB_HAO Maximum AOB growth rate 0.78 d− 1 0.39 / 1.17 qAOB_AMO Maximum rate for the AMO reaction 5.2 g N (g COD) − 1 d1 2.6 / 7.8 KO2_AOB1 AOB affinity constant for O2 (AMO reaction) 1 g O2 m− 3 0.5 / 1.5 KNH4_AOB AOB affinity constant for NH4+ 0.2 g N m− 3 0.1 / 0.3 KO2_AOB2 AOB affinity constant for O2 (HAO reaction) 0.6 g O2 m− 3 0.3 / 0.9 KNH2OH_AOB AOB affinity constant for NH2OH 0.3 g N m− 3 0.15 / 0.45 qAOB_HAO Maximum rate for HAO reaction 5.2 g N (g COD) − 1 d− 1 2.6 / 7.8 KNO_AOB_HAO AOB affinity constant for NO (from HAO) 0.0003 g N m− 3 0.00015 / 0.00045 qAOB_N2O_NN Maximum N2O production rate by NN pathway 0.0078 g N (g COD) − 1 d− 1 0.004 / 0.012 KNO_AOB_NN AOB affinity constant for NO (from NirK) 0.008 g N m− 3 0.004 / 0.012 KO2_AOB_ND AOB constant for O2 effect on the ND pathway 0.5 g O2 m− 3 0.25 / 0.75 KI_O2_AOB N2O constant for production inhibition by O2 0.8 g O2 m− 3 0.4 / 1.2 KHNO2_AOB AOB affinity constant for HNO2 0.004 g N m− 3 0.002 / 0.006 qAOB_N2O_ND Maximum N2O production rate by the ND pathway 1.3 g N (g COD) − 1 d− 1 0.65 / 1.95 KALK_AOB Saturation coefficient for alkalinity (HCO3–) 0.1 mole HCO3– m− 3 0.05 / 0.15 KP_AOB Saturation coefficient for PO43− (nutrient) 0.01 g P m− 3 0.005 / 0.015 µNOB Maximum NOB growth rate 0.78 d− 1 0.39 / 1.17 KO2_NOB Half-saturation coefficient for O2 1.2 g O2 m− 3 0.6 / 1.8 KALK_NOB Saturation coefficient for alkalinity (HCO3–) 0.1 mole HCO3– m− 3 0.05 / 0.15 KNO2_NOB Saturation coefficient for NO2– 0.5 g N m− 3 0.25 / 0.75 KP_NOB Saturation coefficient for PO43− (nutrient) 0.01 g P m− 3 0.005 / 0.015 bAOB Decay rate of AOB 0.096 d− 1 0.048 / 0.144 bNOB Decay rate of NOB 0.096 d− 1 0.048 / 0.144 2.4. ASM2d-N2O structure of all parameters involved in ASM2d-N2O, with calibrated parameters, stoichiometric matrix and kinetics rates can be found in section S9 of the The ASM2d-N2O kinetic model [16] was calibrated to describe the SI. In addition, proportional controllers were modelled in each aerobic Girona WWTP. The ASM2d-N2O model is able to predict COD, N and P reactor to maintain the DO concentration at the measured DO in the removal and N2O production. The model structure is based on the WWTP. The manipulated variable of each proportional controller was ASM2d model developed by Henze et al. [10], and extended to account the mass transfer coefficient for oxygen (kLaO2) in each aerobic reactor. for N2O production with the 2–pathway model for N2O emissions by AOB, developed by Pocquet et al. [14] and the denitrification processes 2.5. Global sensitivity analysis (GSA) with the activated sludge model for nitrogen (ASMN) developed by Hiatt and Grady [27]. The temperature dependence of the biological reactions A GSA was performed to identify the input factors (i.e. parameters) was implemented following the guidelines of the activated sludge model that most affected the model outputs and, therefore, the parameters that No. 2 (ASM2) [10] to describe the different seasonal patterns. The subset should be calibrated preferentially. The selected model output was the 4 B. Solís et al. Chemical Engineering Journal 435 (2022) 134733 calibration cost function (CCF), which is the sum of the squared differ­ reactors. During the preliminary calibration, the influent and opera­ ences between experimental data and dynamic model output (see sec­ tional dynamic data from the experimental campaign were averaged and tion 2.6). used as model inputs (constants inputs). The phosphate concentration in Among the different GSA methods, the Monte Carlo (MC) filtering or the reactors were also averaged and used as output variables. The so­ regional sensitivity analysis (RSA) was a suitable method to select the dium aluminate addition to the IR stream was calibrated using the parameters that were not only more sensitive to CCF but also reduced it XMeOH state variable of ASM2d-N2O, which stands for ferric-hydroxide. [28,29]. RSA is based on mapping the input factors space according to Equation (4) was used as the preliminary calibration cost function whether the associated output, i.e. the CCF, is below (i.e. “behavioural” (PCCF): samples) or above (i.e. “non-behavioural” samples) a predefined √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ √ 7 threshold [30,31]. The workflow used to apply the RSA method was √∑ ( )2 PCCF = √ yexp,i − ymodel,i (4) [29]: 1) a range was defined for the input factor space and a MC i=1 experiment was performed. 2) The model outputs were classified as behavioural (B) or non-behavioural (B) according to the specified Where: “i” is related to each sample in the biological reactor, yexp is the threshold of the CCF and associated to the input factors values. 3) A set averaged experimental phosphate concentration and ymodel is the steady of binary elements was defined, distinguishing between two subsets for state phosphate concentration obtained after a simulation of 300 days. each parameter (Xi): the behavioural subset (Xi |B) and the non- Equation (5) was used to calculate the individual calibration cost behavioural subset (Xi |B). 4) The Smirnov test (equation (3)) was per­ function defined for each monitored variable (CCFi) and Equation (6) formed for each input factor and used as a measure of the Sensitivity was used to calculate the overall CCF: index (Si) [29]; The parameters were ranked in order of influential on √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ CCF reduction by Si. ∑7 √∑ √ n ( )2 CCFi = √ yexpi,j − ymodeli,j (5) Si = max|F(Xi |B) − F(Xi |B) | (3) r=1 j=1 Where F(Xi |B) and F(Xi |B) are the empirical cumulative distribution ∑ 4 functions of the parameter (Xi) when considering the input samples CCF = CCFi = CCFNH4 + CCFNO2 + CCFNO3 + CCFN2O (6) associated with the behavioural and non-behavioural outputs, i=1 respectively. Where: “i” is related to the output variable of interest (NH4+, NO2–, NO3– A total of 59 parameters were included in the GSA study. The or N2O); “j” is related to each experimental data point (n measures); “r” included parameters and their uncertainty ranges are shown in Table 1. is related to each sample zone and yexp and ymodel are related to exper­ The parameters were assumed to be uniformly distributed and the un­ imental data and model output, respectively. The phosphorus related certainty ranges were set according as proposed by Brun et al. [32]. All variables were not added in the CCF for the reasons discussed in section parameters included in the GSA study were kinetic parameters of 2.5. ASM2d-N2O. The hydraulic WWTP parameters were not considered Each dynamic simulation started with a 300-day steady-state (SS) because they were calculated during the hydraulic characterization. The simulation (constant inputs). Then, a 3-day dynamic simulation was influent characterization parameters were neither included, as they performed using the SS simulation results as the initial point and the CCF were measured during the experimental data campaign. Moreover, the was calculated using the operational data sampled during the experi­ stoichiometric parameters of ASM2d–N2O were assumed to be accu­ mental campaign. Each parameter subset was calibrated by minimising rately known parameters and were not included in the GSA. Finally, as the CCF through a global searching minimization method using the the Girona WWTP removes P by chemical precipitation with sodium Matlab function patternsearch. The identifiability of the subsets was aluminate dosage, the parameters related to Polyphosphate Accumu­ evaluated using the FIM methodology [33,34], which was applied to lating Organisms (PAO) were not included. calculate the confidence interval for each parameter, D criteria (CritD), modified E criteria (CritModE), normalized D criteria (NormD) and the 2.6. Calibration procedure ratio between NormD and CritModE (RDE), as detailed in Section S6 of the SI. The ASM2d-N2O model was calibrated once the hydraulics of the The goodness of the model fit for N2O emissions was quantitatively Girona WWTP had been identified and the influent of the experimental assessed with the performance metrics detailed in equations (7) to (10) data campaign had been characterized. The overall procedure followed [35], where yi , y, ̂ y i and n refer to the experimental values, their to calibrate the model is summarized below: average, the model prediction and the number of measurements, 1) a preliminary calibration was performed, aiming to fit the P respectively. chemical removal by sodium aluminate addition. This calibration was ∑n performed under pseudo-steady state conditions to decrease the (̂y i − y)2 R2 = ∑i=1n 2 (7) computational cost. Thus, the experimental values collected during the i=1 (yi − y) experimental campaign were averaged. 2) The CCF was built with the √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅ ∑n dynamic data of nutrients and GHG emissions collected during the y i )2 i=1 (yi − ̂ experimental data campaign. 3) The GSA was performed with the kinetic RMSE = (8) n parameters of the ASM2d-N2O model as input factors and the CCF as output. 4) The top five ranked parameters of the GSA were selected to 1∑n MAE = yi| |yi − ̂ (9) build different subsets of parameters. To minimise the potential problem n i=1 of model overfitting, the parameter subset size to be calibrated was set to a maximum of four parameters. 5) All the subsets were calibrated under 1∑n Bias = yi) (yi − ̂ (10) dynamic conditions. The subset of the optimized parameters that most n i=1 reduced the CCF was selected. The identifiability of each subset and the The experimental campaign lasted three days. In addition, the gas confidence interval for its parameters were evaluated based on the FIM collection system could only analyse one aerobic zone at a time and, methodology. unfortunately, the system failed during one night due to the gas mea­ The aim of the preliminary calibration was to fit the sodium alumi­ surement computer was powered off. Therefore, in the absence of a large nate addition to describe the phosphate concentration in the biological 5 B. Solís et al. Chemical Engineering Journal 435 (2022) 134733 dynamic data set, all available data from this experimental campaign were used for the calibration process. The validation was performed with data selected from a previous N2O measuring campaign in the same WWTP [21]. Specifically, the data used for model validation were the N2O emissions and the ammonium profiles measured in the AER2 compartment of the WWTP during 15–21 November 2016 and reported in a previous experimental work [21]. The validation procedure and results are detailed in section S8 of the SI. 3. Results and discussion 3.1. Experimental data and plant performance The NH4+, NO2–, NO3– and PO43− profiles in the influent and in each compartment are shown in Figures S2 to S5 in section S2 of the SI. The influent temperature and pH were approximately constant during the experimental campaign at 24.4 ± 0.5 ◦ C and 7.44 ± 0.09, respectively. Good plant performance was achieved, obtaining a COD removal of 96% and higher than 99% for total Kjeldahl nitrogen (TKN) and P. Mean DO values were 1.8, 1.5, 1.2 and 2.0 g O2 m− 3 in AER1 to AER4 compart­ ments, respectively. Nitrite concentrations were below 0.25 g N m− 3 in all compartments, with averaged ammonium and nitrate concentrations in the last compartment of 0.1 g N m− 3 and 3.3 g N m− 3, respectively, showing high AOB and NOB activity. Figure S5 shows that the PO43− concentration only increased on average by 2.5 g P m− 3 in the anaerobic reactor (ANA1), compared to the PO43− concentration in the influent, showing a low PAO activity. Hence, the addition of chemical P precip­ itant probably limited the PAO activity. The COD fractionation characterization (i.e. the calculation of the COD model state variables from the measured variables) was performed on the basis of the effluent filtered COD (CODfil,eff) and three different influent COD measurements: total COD (CODtot), filtered COD at 1.2 µm (CODfil) and flocculated-filtered COD at 0.45 µm (CODsol). More details on this procedure, the list of all determined model influent state vari­ ables and the COD influent profiles can be found in section S3 of the SI. The influent fractionation was assumed constant during the dynamic calibration. Fig. 2. Measured N2O emissions and NH4+ concentration in compartments Fig. 2 shows the N2O emissions from the first three aerated zones AER1, AER2 and AER3. The grey area represents a technical failure in the N2O (AER1 to AER3) together with the ammonium concentration profile gas measurement system. The dashed line at the N2O and NH4+ analytical data obtained from analytical tests and from an online sensor in AER2. The points is an aid for better visualization of the experimental profiles. grey area in Fig. 2 indicates that N2O data was not available due to a technical failure. Similar peak profiles of N2O emissions were observed 3.2. Hydraulic characterization for the other three aerobic reactors monitored. However, the amount of N2O emitted was different, with AER1 being the compartment with The pulse response at the influent of biological reactors and the highest emissions and AER3 with the lowest. This could be related to the flowrates measured during the RTD experiment are shown in Fig. 3. The fact that higher ammonium concentrations were present in AER1, that influent tracer concentration, the flowrates and the dimensional data of led to higher ammonium oxidation rates. Furthermore, Fig. 2 shows that the Girona WWTP were used as model inputs to characterize the hy­ N2O emissions were related to ammonium concentration and thus draulics of the plant. emissions decreased to negligible levels when ammonium was depleted. Fig. 4 shows the experimental bromide concentrations measured The same pattern of N2O peak emissions were reported in other studies during the tracer test. The concentration of Br− before the KBr pulse was in full-scale WWTPs or in lab experiments with nitrifying-enriched not negligible. Therefore, a constant inlet bromide concentration of sludge [26,36]. The peak is attributed to the sudden increase of 0.0627 g m− 3 was considered during the experiment. A tracer mass ammonium, which produces a transient between low-activity to high- balance over the secondary treatment after 78 h showed that a 79% of activity of nitrifying biomass [4,9,26]. During the experimental the total tracer introduced was detected at the output and that 3.5 kg of campaign, the averaged N2O emission factor (N2O-EF) (calculated as the the injected Br− remained in the reactor when the experiment was percentage of the influent TKN load emitted as N2O-N), of AER1 to AER3 stopped. reactors was 0.41%, which is in the low range of the N2O-EF reported for The initial assumption that each reactor operated as an ideal CSTR full-scale WWTPs [4,9,12,26]. However, the calculated N2O-EF was was correct, as the trends of the model predictions agreed with the slightly higher than the N2O-EF calculated during a large monitoring experimental values. However, the assumption that the secondary campaign at the same WWTP, N2O-EF of 0–0.13% [26]. The increase in settler flow pattern was also a CSTR was false and therefore the model N2O-EF can be due to the multiple factors that can be found in an urban was revised. The parameters to be estimated were the number of N WWTP, such as operational parameters (e.g. temperature, DO, pH and tanks-in-series of the secondary settler and the percentage distribution recycling flows) and changes in wastewater flow and composition (e.g. of the influent flow between each biological line (i.e. fQ1 and fQ2, where nitrogen concentration, COD concentration and fractionation, and in­ fQ1 and fQ2 are the percentage of the influent flow going to the first and crease in the fraction of industrial wastewater) [12]. second biological treatment lines, respectively, and fQ1 + fQ2 = 1). 6 B. Solís et al. Chemical Engineering Journal 435 (2022) 134733 Fig. 3. Pulse response concentration measured in the effluent of the primary clarifier (influent of the biological reactors) and flowrates measured during the tracer experiment. Fig. 4. Tracer experiment results and model results in each measured compartment and in the effluent. Fig. 4 shows the results obtained during the hydraulic calibration. 3.3. Preliminary calibration The model accurately described the trends of the experimental data. The optimized parameters for the secondary settler flow pattern were N = 5 The next step after the identification of the WWTP hydraulics was the tanks and hence the secondary settler flow pattern was closer to a plug calibration of the kinetic model ASM2d-N2O. As the flow distribution flow reactor flux model than to a CSTR. The flowrate was distributed between both biological treatment lines was approximately the same, approximately equally between the two lines, with 49.33% of the only one biological treatment line was considered in the kinetic cali­ influent flowrate going to the first biological treatment line and 50.67% bration and in the GSA studies. The preliminary calibration (see section to the second treatment line. These results contradicted the initial 2.6) aimed to reduce structural discrepancies between the model and the thinking of the plant operators that the distribution was not experimental variables [37], particularly for those related to P, before symmetrical. performing the GSA and the subsequent dynamic calibration. 7 B. Solís et al. Chemical Engineering Journal 435 (2022) 134733 The XMeOH addition that minimized the PCCF was 184 kg Fe(OH)3 Table 2 d− 1. The PCCF obtained was 2.7, achieving a reduction of 74% Ranking of the first twenty parameters obtained in the GSA. compared to the PCCF before the preliminary calibration. The fits ob­ Position Parameter Description tained on the average phosphate concentration are shown in Figure S7 of 1 µNOB Maximum NOB growth rate section S4 of the SI. With this addition of XMeOH, the model was found to 2 bNOB Decay rate of NOB predict that the percentage of P biologically removed by PAO was 24%. 3 KNO2_NOB Saturation coefficient for NO2– The calibrated XMeOH addition value was maintained constant during all 4 qAOB_AMO Maximum rate for the AMO reaction dynamic simulations. Moreover, as phosphate concentrations were 5 KO2_NOB Half-saturation coefficient for O2 6 nG5 Anoxic growth factor (N2O → N2) calibrated by chemical precipitation, PAO-related parameters were not 7 KOH5 Half-saturation coefficient for O2 considered during the ASM2d–N2O calibration and phosphate concen­ 8 qAOB_HAO Maximum rate for HAO reaction trations were not considered in the CCF, i.e. the PAO-related parameters 9 qAOB_N2O_ND Maximum N2O production rate by the ND pathway were maintained at their default values in all dynamic simulations. 10 KNO2 Saturation/inhibition coefficient for NO2– 11 KI5NO NO inhibition coefficient (N2O → N2) 12 KI_O2_AOB N2O constant for production inhibition by O2 3.4. Development of the calibration cost function (CCF) 13 KO2_AOB1 AOB affinity constant for O2 (AMO reaction) 14 nG3 Anoxic growth factor (NO2–→NO) The data obtained during the experimental campaign were used to 15 KS5 Half-saturation coefficient for substrate dynamically calibrate the ASM2d-N2O model to the Girona WWTP. In 16 KNO2_Den Half-saturation coefficient for NO2– 17 KI4NO NO inhibition coefficient (NO → N2O) total, four variables were included in the cost function and considered as 18 KALK Saturation coefficient for alkalinity (HCO3–) the output variables of interest: ammonium, nitrite and nitrate in the 19 KHNO2_AOB AOB affinity constant for HNO2 biological reactors and N2O emissions from AER1 to AER3. Therefore, 20 KNO3_H Saturation/inhibition coefficient for NO3– the CCF was divided into four different cost functions (CCFi) for each output variable of interest (i.e. CCFNH4, CCFNO2, CCFNO3, CCFN2O). The 3.6. Dynamic calibration overall CCF (Equation (6)) was calculated as the sum of each individual output variable CCFi (Equation (5)). As all the CCFi were of the same The dynamic calibration was conducted after identifying the ASM2d- order of magnitude, no weighting factors were included in the CCF N2O parameters most likely to reduce the CCF. Different parameter calculation. subsets were defined with all the possible combinations of the five top ranked parameters from the GSA (Table 2). The size of parameter subsets 3.5. Global sensitivity analysis ranged from one to four parameters, resulting in a total of 30 parameter subsets to be calibrated (Table S3 in SI). The parameter subset and A number of N = 2000 simulations was selected for the MC experi­ calibrated values that most reduced the CCF were µNOB = 0.674 ± ment. The CCF and each CCFi for each simulation were then evaluated 0.019 d–1, qAOB_AMO = 5.517 ± 0.058 gN (gCOD)–1 d–1, KO2_NOB = 0.126 and discretized into two populations, B and B (see section 2.5). The ± 0.029 gO2 m− 3 and KNO2_NOB = 0.126 ± 0.030 gN m− 3. The CCF was threshold fixed to discretize the CCF (and the associated input param­ reduced by 53.3% with this parameter subset, compared to the CCF after eters) was the maximum reduction on the CCF that can be achieved with preliminary calibration, mainly due to the reduction of the CCF of NO2– a number of simulations of the behavioural group of at least NB = 100, and N2O (87.7 and 86.5%, respectively). The fit between experimental representing 5% of all simulations. The maximum reduction in CCF that data and model predictions for N-species is shown in Fig. 5 and those for satisfied the requirement of NB = 100, and therefore, the threshold N2O emissions are shown in Fig. 6. selected for GSA evaluation was 40%, compared to the CCF calculated The calibration results for each subset of parameters can be found in with the default parameter values of ASM2d-N2O. For this selected Table S3 of the SI. The selected subset shows low confidence interval threshold, the top ranked GSA indices are shown in Table 2. Table S2 in values for all parameters and a parameter confidence interval norm of section S5 of the SI shows the full rank of parameters evaluated. 33.6%, which is the lowest of all subsets with four parameters. It also Among the CCFi-related results, a high reduction threshold of around provides a high NormD = 9.8 1017, a low CritModE = 809 and the 80% was found for nitrite and nitrous oxide CCFs. On the other hand, the highest RDE = 1.22⋅1015. These values, compared to those of the other thresholds found for the ammonium and nitrate CCFs were around 10%. subsets, indicate that the selected subset has a good identifiability linked These results showed that the maximum reduction in the overall CCF to the ability to provide a low CCF. In contrast, the subset µNOB, bNOB, that could be achieved was related to the reduction in the nitrite and qAOB_AMO, KO2_NOB also provides a very similar CCF, but it has an nitrous oxide CCFs. The initial simulations agreed with this observation, extremely high parameter confidence interval for bNOB (1211%), i.e. as the predicted ammonium and nitrate concentrations trends agreed bNOB = 0.0008 ± 0.0095 d–1, which shows the low identifiability of this with the experimental values. In addition, Table 2 shows that the five parameter. The problem of a much higher confidence interval than the top ranked parameters were related to nitrification processes, more parameter value with bNOB also appears in all subsets where very low specifically to the nitratation process, as four of the five top ranked bNOB values are obtained during the optimization. This indicates that, parameters were related to NOB (µNOB, bNOB, KNO2_NOB and KO2_NOB) and with the experimental information available in this study, it would not the other parameter was related to AOB bacteria (qAOB_AMO). Therefore, be advisable to use bNOB for model fitting. In fact, this is the only the GSA results showed that the ASM2d-N2O model calibration for the parameter of the initial five that is not included in the chosen subset. Girona WWTP should focus on the calibration of the nitrifying bacteria Regarding the predictions obtained, the model reasonably predicts parameters, which is not surprising. However, in view of good modelling the concentration trends of the different N-species in the reactor. The practices, a systematic methodology to understand the parameter ASM2d-N2O was able to explain the low nitrite concentrations measured sensitivity should always be incorporated, particularly when there are a during the experimental campaign, by reducing the NOB oxygen and high number of parameters. In this case, i.e. for this model, specific nitrite affinity constants with respect to the default values [16]. Exper­ WWTP operational conditions and available experimental data, these imental nitrate data and model fits showed a discrepancy in nitrification results allowed to detect the most sensitive parameter among all those capacity, as measured nitrate showed smaller increases among aerobic related to the biological nitrogen removal process and confirmed that reactors (AER1 to AER3) compared to model predictions. The model any unexpected parameter had a high sensitivity on N2O measurements. could not predict such a high degree of simultaneous nitrification and denitrification occurring in AER2 and AER3. This was due to the low DO levels in AER2 and AER3 (1.5 and 1.2 g O2 m− 3, respectively) in addition 8 B. Solís et al. Chemical Engineering Journal 435 (2022) 134733 Fig. 5. Concentrations of ammonium, nitrite and nitrate measured during the experimental campaign and fit obtained during the dynamic calibration of ASM2d-N2O. to the inherent instability of the DO control system in the plant, which concentration in AER1, obtaining an average value of 0.17 g N m− 3 caused significant noise in the DO concentration and favoured denitri­ during the three days of dynamic simulation. This a priori low nitrite fication processes when the DO concentration was low. concentration was high enough to support the ND pathway in this On the other hand, the goodness of the ASM2d-N2O fit for N2O reactor, thus causing the overestimation of N2O emissions. experimental data was evaluated in more detail. Statistical model fit ND is the biological N2O production pathway responsible for most of criteria were calculated for each reactor (Table 3) as detailed in section emissions during wastewater treatment [9,12,16,38]. In fact, the 2.6. The coefficient of determination was low in the three reactors (R2 = average contribution of the ND pathway to the total N2O production in 0.107, 0.518 and 0.533), indicating that the model could only partially AER1 was 82% and decreased to 78% and 48% in AER2 and AER3, describe the experimental results. The highest RMSE value resulted in respectively (see Figure S8 in SI). These values are in agreement with AER1, but when comparing the RMSE with the average experimental literature ranges for aerobic reactors [38]. The NN pathway contributed value (y), the error was 85%, 79% and 139% for reactors AER1, AER2 6% to the total N2O production and the HD pathway contributed mostly and AER3 respectively. The MAE results were similar to those of RMSE, in AER3 with 45% (Figure S8 in SI), where the DO was the lowest of the obtaining the highest value for AER1. Finally, the model bias was aerobic zones (1.2 g O2 m− 3). The higher contribution of the HD negative for all three reactors (Bias = 0.114, -0.285 and -0.120), indi­ pathway in AER3 shows that simultaneous nitrification and denitrifi­ cating a slight model underestimation of N2O emissions. Despite the cation occurred in this reactor due to the low DO coupled with low NH4+ statistical values discussed above, the model was able to reasonably and high NO3– concentrations. The predicted N2O-EF, only considering capture the dynamics and range of N2O emissions in AER2 and AER3 AER1 to AER3 zones, was 0.55% which is very similar to the measured (Fig. 6), except for the period 1.0–1.5 d, where unexpectedly low N2O one (0.41%). values were registered by the online measuring system in all three re­ Although the model predictions are in line with previous work in the actors. These reasonable predictions were obtained even considering literature, it should be noted that biological N2O production are complex that only one parameter related to biomass producing N2O was modified processes with several metabolic pathways with challenging in­ (qAOB_AMO) and also taking into account the intrinsic operational vari­ teractions that have probably not been fully revealed. There is a lack of a ability observed in full-scale WWTPs. consensus model describing all the scenarios and, thus, the biological Going into more detail with respect to the discrepancy of N2O model used may not be capable to describe all the particular cases. The emissions in the period 1.0–1.5d, the model predicted a peak of ASM2d-N2O model contains several assumptions and short-cuts in order ammonium during this period that was not experimentally detected to describe full-scale N2O emissions in a simple way. For instance, the (Fig. 5). The model predictions of N2O emissions during this second ASM2d-N2O model assumes that the AOB-driven nitrite reduction to daily peak of ammonium were mostly related to the ND pathway, since, N2O is a one-step process to avoid a NO loop [16] based on the original as a result of the ammonium accumulation, the model also predicted a structure [14] that was only focused on AOB metabolism. Likewise, the slight accumulation of hydroxylamine up to 0.9 g N m− 3 (three times the quantity and quality dataset used for calibration corresponds to a rela­ affinity constant of AOB), which is the electron acceptor substrate for the tively short-term period and, as such, we missed rain-events or periods ND pathway. Moreover, the model predicted the highest nitrite with changing temperature or variable flow. Therefore, these 9 B. Solís et al. Chemical Engineering Journal 435 (2022) 134733 distribution was simulated, i.e. 50% of the influent flowrate was fed to each line, according to the tracer experiments results (section 3.2). The second case was modelled assuming that 40% of the influent flowrate went to the first treatment line and the remaining 60% to the second one. This study was performed to assess the possible effect than an incorrect flow distribution between treatment lines had on the N2O emissions due to possible treatment line overloading and changes in treatment plant loads. Each simulation was done following the same methodology as in the dynamic calibration (section 2.6) and with the same model inputs. In addition, it was assumed that the aeration system was able to maintain the same DO concentration despite the increase of N load. Fig. 7 shows the predicted N2O emissions in each treatment line and the total N2O emissions and N2O-EF for each case. Fig. 7A shows that the predicted N2O emissions from the treatment lines increased with increasing influent flowrate. The average N2O emission rate obtained for an influent flowrate distribution of 40, 50 and 60% were 2.5, 3.2 and 3.8 kg N2O-N d− 1 respectively. The same ammonium was obtained in the effluent for all three influent distribution simulations, however, nitrite and nitrate in the effluent increased with increasing influent flowrate due to the increase on the TKN influent load. Fig. 7A reveals that the same total N2O emissions were predicted for both cases (6.3 kg N2O-N d− 1). This is because the N2O emissions increased linearly with increasing influent flowrate and therefore the total N2O emissions for both cases were the same. Fig. 7B also shows that the predicted N2O-EF for both lines was the same for each case because the N2O emissions increased with the same slope as the TKN removed. Therefore, these simulations show that if one of the two lines was overloaded, emissions would increase in that line, but these higher emissions would be compensated to the same extent by the lower emissions in the under­ loaded line. Thus, the plant, would not suffer an increase in total N2O emissions even if the flow distribution was not entirely correct, as long Fig. 6. N2O emissions measured during the experimental campaign and fits obtained during the dynamic calibration of ASM2d-N2O. In the grey area no experimental N2O emission data were available due to a technical failure. Table 3 Performance criteria to assess the goodness of fit for N2O predictions. n R2 y RMSE MAE Bias AER1 reactor 37 0.107 1.339 1.141 1.014 − 0.114 AER2 reactor 37 0.518 0.763 0.603 0.514 − 0.285 AER3 reactor 38 0.533 0.289 0.401 0.310 − 0.120 considerations, i.e. the model structure and the quantity/quality of the calibration dataset, may have a strong influence on the parameter subset selection, parameter estimates and model predictions, and should be taken into account when analysing the results obtained, considering also the possible model limitations when it is to be applied to other scenarios. 4. Model exploitation Once the dynamic calibration and the model validation (Section S8 Fig. 7. Predicted N2O emissions in two different flowrate scenarios: Case 1 in SI) were done, two scenarios were simulated to investigate the effect (equal influent flowrate distribution in both treatment lines) and Case 2 on N2O emissions and N2O-EF of varying the flowrate distribution be­ (40–60% influent flowrate distribution). (A) N2O emissions in each treatment line for 40%-50%-60% flowrate and total N2O emissions. (B) Predicted N2O-EF tween both treatment lines (Fig. 1). In the first case, an equal flow for both cases. 10 B. Solís et al. Chemical Engineering Journal 435 (2022) 134733 as it can be ensured that the aeration system is capable of providing References sufficient oxygen to maintain the DO around the setpoint. For scenarios where there is no oxygen control and constant aeration is used, an in­ [1] IPCC, 2014: Climate Change 2014: Synthesis Report . Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel crease in ammonium loading may decrease the DO concentration, on Climate Change [Core Writing Team, R.K. Pachauri and L.A. Meyer (eds.)]. moving the operating point to one with higher N2O production, which IPCC, Geneva, Switze, (n.d.). would increase emissions as previously studied [16]. [2] D.J.I. Gustavsson, S. Tumlin, Carbon footprints of Scandinavian wastewater treatment plants, Water Sci. Technol. 68 (2013) 887–893, https://doi.org/ 10.2166/wst.2013.318. 5. Conclusions [3] A.R. Ravishankara, J.S. Daniel, R.W. Portmann, Nitrous Oxide (N 2 O): The Dominant Ozone-Depleting Substance Emitted in the 21st Century, Science 326 (5949) (2009) 123–125, https://doi.org/10.1126/science.1176985. This work is a comprehensive calibration of the ASM2d-N2O model [4] M.J. Kampschreur, H. Temmink, R. Kleerebezem, M.S.M. Jetten, M.C.M. van to a full-scale WWTP including hydraulics with the following main Loosdrecht, Nitrous oxide emission during wastewater treatment, Water Res. 43 findings: (17) (2009) 4093–4103, https://doi.org/10.1016/j.watres.2009.03.001. [5] J. Foley, D. de Haas, Z. Yuan, P. Lant, Nitrous oxide generation in full-scale biological nutrient removal wastewater treatment plants, Water Res. 44 (3) (2010) - Modelling the flow patterns in the plant is very interesting in view of 831–844, https://doi.org/10.1016/j.watres.2009.10.033. its calibration. The tracer experiment showed that all reactors of the [6] B.-J. Ni, Z. Yuan, Recent advances in mathematical modeling of nitrous oxides two treatment lines had a correct hydraulic behaviour, as no dead emissions from wastewater treatment processes, Water Res. 87 (2015) 336–346, https://doi.org/10.1016/j.watres.2015.09.049. volumes, flux recycling or by-passes were found. Furthermore, it was [7] J.D. Caranto, K.M. Lancaster, Nitric oxide is an obligate bacterial nitrification demonstrated that an equal flow was flowing to each line. intermediate produced by hydroxylamine oxidoreductase, Proc. Natl. Acad. Sci. - A GSA method (RSA) was successfully applied to rank the parameters 114 (31) (2017) 8217–8222, https://doi.org/10.1073/pnas.1704504114. [8] P. Wunderlin, J. Mohn, A. Joss, L. Emmenegger, H. Siegrist, Mechanisms of N2O most likely to reduce the CCF. The sensitivity analysis revealed that production in biological wastewater treatment under nitrifying and denitrifying the top ranked parameters were related to nitrifying organisms. conditions, Water Res. 46 (4) (2012) 1027–1037, https://doi.org/10.1016/j. - The dynamic kinetic calibration of ASM2d-N2O correctly captures watres.2011.11.080. [9] Y. Law, L. Ye, Y. Pan, Z. Yuan, Nitrous oxide emissions from wastewater treatment the trends of N-species in the reactors, only by modifying four kinetic processes, Philos. Trans. R. Soc. B Biol. Sci. 367 (1593) (2012) 1265–1277, https:// parameters. doi.org/10.1098/rstb.2011.0317. - The N2O-EF predicted was very similar to that measured experi­ [10] M. Henze, W. Gujer, T. Mino, M.C.M. van Loosdrecht, Activated Sludge Models ASM1, ASM2, ASM2d and ASM3, Scientific and technical report No 9, IWA mentally and the emission profiles trends were in general in agree­ Publishing (2000), https://doi.org/10.2166/9781780402369. ment with the experimental data. [11] B.-J. Ni, L. Ye, Y. Law, C. Byers, Z. Yuan, Mathematical modeling of nitrous oxide - The simulations revealed that, under current plant operating condi­ (N2O) emissions from full-scale wastewater treatment plants, Environ. Sci. Technol. 47 (14) (2013) 7795–7803, https://doi.org/10.1021/es4005398. tions, an increase in flow in one of the treatment lines would increase [12] T.M. Massara, S. Malamis, A. Guisasola, J.A. Baeza, C. Noutsopoulos, E. Katsou, N2O emissions in that overloaded line, but would be offset by a A review on nitrous oxide (N2O) emissions during biological nutrient removal from decrease in the underloaded line, so that total emissions would not municipal wastewater and sludge reject water, Sci. Total Environ. 596-597 (2017) vary significantly, provided the control system was able to maintain 106–123, https://doi.org/10.1016/j.scitotenv.2017.03.191. [13] B.-J. Ni, Z. Yuan, K. Chandran, P.A. Vanrolleghem, S. Murthy, Evaluating four the operating DO. mathematical models for nitrous oxide production by autotrophic ammonia- oxidizing bacteria, Biotechnol. Bioeng. 110 (1) (2013) 153–163, https://doi.org/ Declaration of Competing Interest 10.1002/bit.24620. [14] M. Pocquet, Z. Wu, I. Queinnec, M. Spérandio, A two pathway model for N2O emissions by ammonium oxidizing bacteria supported by the NO/N2O variation, The authors declare that they have no known competing financial Water Res. 88 (2016) 948–959, https://doi.org/10.1016/j.watres.2015.11.029. interests or personal relationships that could have appeared to influence [15] C. Domingo-Félez, B.F. Smets, A consilience model to describe N2O production during biological N removal, Environ. Sci. Water Res. Technol. 2 (6) (2016) the work reported in this paper. 923–930, https://doi.org/10.1039/C6EW00179C. [16] T.M. Massara, B. Solís, A. Guisasola, E. Katsou, J.A. Baeza, Development of an Acknowledgements ASM2d-N2O model to describe nitrous oxide emissions in municipal WWTPs under dynamic conditions, Chem. Eng. J. 335 (2018) 185–196, https://doi.org/10.1016/ j.cej.2017.10.119. Borja Solís is grateful for the PIF PhD grant funded by Universitat [17] M. Gresch, D. Braun, W. Gujer, The role of the flow pattern in wastewater aeration Autònoma de Barcelona. Borja Solís, Albert Guisasola and Juan Antonio tanks, Water Sci. Technol. 61 (2010) 407–414, https://doi.org/10.2166/ wst.2010.803. Baeza are members of the GENOCOV research group (Grup de Recerca [18] D. Olivet, J. Valls, M.À. Gordillo, À. Freixó, A. Sánchez, Application of residence Consolidat de la Generalitat de Catalunya, 2017 SGR 1175, www. time distribution technique to the study of the hydrodynamic behaviour of a full- genocov.com) and were supported by the European Union’s Horizon scale wastewater treatment plant plug-flow bioreactor, J. Chem. Technol. 2020 research and innovation program under the Marie Skłodowska- Biotechnol. 80 (4) (2005) 425–432, https://doi.org/10.1002/jctb.1201. [19] F. Wang, J. Li, D. Bian, Z. Nie, S. Ai, X.i. Tian, H. Wang, Effect of hydraulic Curie C-FOOT-CTRL project (grant agreement No 645769). M. Pijuan residence time and inlet flow distribution ratio on the pollutant removal of low- and Ll. Corominas acknowledge the support from the Economy and temperature municipal wastewater in multistage AO process, Water Environ. J. 34 Knowledge Department of the Catalan Government through a Consoli­ (S1) (2020) 742–752, https://doi.org/10.1111/wej.12574. [20] F. Sánchez, A. Viedma, A.S. Kaiser, Hydraulic characterization of an activated dated Research Group (ICRA-TECH - 2017 SGR 1318) - Catalan Institute sludge reactor with recycling system by tracer experiment and analytical models, for Water Research. A special thanks to Ms. Cristina Pujol and Mr. Lluís Water Res. 101 (2016) 382–392, https://doi.org/10.1016/j.watres.2016.05.094. Ayach from TRARGISA, the company responsible for the operation of [21] M.C. Collivignarelli, M. Carnevale Miino, S. Manenti, S. Todeschini, E. Sperone, G. Cavallo, A. Abbà, Identification and localization of hydrodynamic anomalies in the WWTP for all the help during the monitoring campaign. We a real wastewater treatment plant by an integrated approach: RTD-CFD Analysis, acknowledge the contribution of Dr. Ayla Kiser (www.aylakiser.com) in Environ. Process. 7 (2) (2020) 563–578, https://doi.org/10.1007/s40710-020- the design and execution of the tracer test. 00437-4. [22] M. Majewsky, T. Gallé, M. Bayerle, R. Goel, K. Fischer, P.A. Vanrolleghem, Xenobiotic removal efficiencies in wastewater treatment plants: Residence time Appendix A. Supplementary data distributions as a guiding principle for sampling strategies, Water Res. 45 (18) (2011) 6152–6162, https://doi.org/10.1016/j.watres.2011.09.005. [23] G. Tchobanoglous, F. Burton, H.D. Stensel, Wastewater engineering: Treatment and Supplementary data to this article can be found online at https://doi. reuse, Am. Water Work. Assoc. J. 95 (2003) 201. org/10.1016/j.cej.2022.134733. [24] F. Coen, B. Petersen, P.A. Vanrolleghem, B. Vanderhaegen, M. Henze, Model-based characterisation of hydraulic, kinetic and influent properties of an industrial WWTP, Water Sci. Technol. 37 (1998) 317–326, https://doi.org/10.1016/S0273- 1223(98)00352-7. [25] P.A. Vanrolleghem, G. Insel, B. Petersen, G. Sin, D. De Pauw, I. Nopens, H. Dovermann, S. Weijers, K. Gernaey, A comprehensive model calibration 11 B. Solís et al. Chemical Engineering Journal 435 (2022) 134733 procedure for activated sludge models, Proc. Water. Environ. Fed. 2003 (9) (2003) Res. 36 (16) (2002) 4113–4127, https://doi.org/10.1016/S0043-1354(02)00104- 210–237, https://doi.org/10.2175/193864703784639615. 5. [26] A. Ribera-Guardia, L. Bosch, L. Corominas, M. Pijuan, Nitrous oxide and methane [33] V.C. Machado, G. Tapia, D. Gabriel, J. Lafuente, J.A. Baeza, Systematic emissions from a plug-flow full-scale bioreactor and assessment of its carbon identifiability study based on the Fisher Information Matrix for reducing the footprint, J. Clean. Prod. 212 (2019) 162–172, https://doi.org/10.1016/j. number of parameters calibration of an activated sludge model, Environ. Model. jclepro.2018.11.286. Softw. 24 (11) (2009) 1274–1284, https://doi.org/10.1016/j.envsoft.2009.05.001. [27] W.C. Hiatt, C.P.L. Grady, An updated process model for carbon oxidation, [34] A. Guisasola, J.A. Baeza, J. Carrera, G. Sin, P.A. Vanrolleghem, J. Lafuente, The nitrification, and denitrification, Water Environ. Res. 80 (11) (2008) 2145–2156, influence of experimental data quality and quantity on parameter estimation https://doi.org/10.2175/106143008X304776. accuracy, Educ. Chem. Eng. 1 (1) (2006) 139–145, https://doi.org/10.1205/ [28] A. Saltelli, M. Ratto, S. Tarantola, F. Campolongo, Sensitivity analysis for chemical ece06016. models, Chem. Rev. 105 (7) (2005) 2811–2828, https://doi.org/10.1021/ [35] D. Wallach, Evaluating crop models. Work. with Dyn. Crop Model. Eval. Anal. cr040659d. Parameterization, Appl, Elsevier, 2006, pp. 11–54. https://www.elsevier. [29] A. Saltelli, M. Ratto, S. Tarantola, F. Campolongo, Sensitivity analysis practice: A com/books/working-with-dynamic-crop-models/wallach/978-0-444-52135-4. guide to scientific models (2004), https://doi.org/10.1016/j.ress.2005.11.014. [36] A. Ribera-Guardia, M. Pijuan, Distinctive NO and N2O emission patterns in [30] A. Saltelli, M. Ratto, T. Andres, F. Campolongo, J. Cariboni, D. Gatelli, M. Saisana, ammonia oxidizing bacteria: Effect of ammonia oxidation rate, DO and pH, Chem. S. Tarantola, Global Sensitivity Analysis: The Primer (2008), https://doi.org/ Eng. J. 321 (2017) 358–365, https://doi.org/10.1016/j.cej.2017.03.122. 10.1111/j.1751-5823.2008.00062_17.x. [37] V.C. Machado, J. Lafuente, J.A. Baeza, Activated sludge model 2d calibration with [31] F. Pianosi, K. Beven, J. Freer, J.W. Hall, J. Rougier, D.B. Stephenson, T. Wagener, full-scale WWTP data: comparing model parameter identifiability with influent and Sensitivity analysis of environmental models: A systematic review with practical operational uncertainty, Bioprocess Biosyst. Eng. 37 (7) (2014) 1271–1287, workflow, Environ. Model. Softw. 79 (2016) 214–232, https://doi.org/10.1016/j. https://doi.org/10.1007/s00449-013-1099-8. envsoft.2016.02.008. [38] G. Tallec, J. Garnier, G. Billen, M. Gousailles, Nitrous oxide emissions from [32] R. Brun, M. Kühni, H. Siegrist, W. Gujer, P. Reichert, Practical identifiability of secondary activated sludge in nitrifying conditions of urban wastewater treatment ASM2d parameters - Systematic selection and tuning of parameter subsets, Water plants: Effect of oxygenation level, Water Res. 40 (15) (2006) 2972–2980, https:// doi.org/10.1016/j.watres.2006.05.037. 12