Microsoft Word - 8bove.docx CHEMICAL ENGINEERING TRANSACTIONS VOL. 50, 2016 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Guest Editors: Katharina Kohse-Höinghaus, Eliseo Ranzi Copyright © 2016, AIDIC Servizi S.r.l., ISBN 978-88-95608-41-9; ISSN 2283-9216 Wastewater Stabilisation Ponds System: Global Sensitivity Analysis on Network Design Maria Paz Ochoa, Vanina Estrada, Patricia M. Hoch* Planta Piloto de Ingeniería Química – PLAPIQUI –CONICET p.hoch@plapiqui.edu.ar In this work, a wastewater treatment network (WWTN) consisting of a system of ponds is modelled rigorously, taking into account dynamic mass balances for the main groups of bacteria, together with different types of organic load, algae biomass, nutrients, etc. To obtain the optimal configuration, we had first formulated a superstructure embedding these rigorous models as a mixed integer non-linear programming (MINLP) problem, with the objective to synthesise and design the WWTN that minimises the total annual costs, subject to environmental regulations. Then, a global sensitivity analysis (GSA) is performed on this kinetic dynamic model for the obtained optimal configuration of three stabilisation ponds (two aerobic ponds in series followed by a facultative one) to determine the most influential parameters of the model considering the whole range of parameters variation, as well as parameter ranking. GSA is implemented using Sobol’s method, a variance based technique. The technique is implemented within gPROMS platform, a differential algebraic equation oriented environment where stochastic simulations are performed. Temporal profiles for the first order, total order and interactional sensitivity indices are obtained for the main differential and algebraic state variables. Numerical results provide useful information about the complex relationships between technological, economic and environmental variables of the processes in the WWTN design optimisation. 1. Introduction In a previous work, a WWTN embedding rigorous models of different types of stabilisation ponds was developed with synthesis and design purposes by Ochoa et al. (2016a). The model was formulated as a MINLP problem to establish optimal configuration system and pond dimensions necessary to conform to environmental discharge regulations. In addition, different scenarios of inlet concentration of organic load were explored to take into account process uncertainty leading to the conclusion that the model was highly influenced by initial organic load. In this paper, another type of uncertainty is considered, related to time-invariant model parameters such as kinetic parameters, whose uncertainty results from measurement errors or the impossibility of modelling exactly the physical behaviour of a system. The values of state variable can be greatly influenced by this type of uncertainty. For this reason, it is important to identify the parameters to which model state variables are most sensitive, which is achieved in general by a sensitivity analysis (SA). Techniques for sensitivity analysis can be classified into local and global. Local techniques are based on Taylor series expansion, requiring linearity and additivity. Global sensitivity analysis (GSA) aims to quantify the relative importance of input variables or factors in determining the value of an output variable. GSA method should be used when the model is nonlinear and various input variables are affected by uncertainties of different orders of magnitude, because they take into account the influence of the whole range of variation and the form of the probability density function in the input. GSA method computes the effect of factor xi while all others xj, j≠i, are varied as well. Temporal profiles of the influence of factors and input variables are obtained when dynamic models are analysed, leading to a great insight on the sources of uncertainty during the time horizon (Saltelli et al., 2008). The wastewater stabilisation pond network model and the GSA methodology were implemented in an equation oriented environment with a differential algebraic equation solver in gPROMS (PSEnterprise, 2014). The implemented GSA strategy is variance-based (Sobol’, 1993) and allows the determination and DOI: 10.3303/CET1650032 Please cite this article as: Ochoa M.P., Estrada V., Hoch P., 2016, Wastewater stabilisation ponds system: global sensitivity analysis on network design, Chemical Engineering Transactions, 50, 187-192 DOI: 10.3303/CET1650032 187 classification of model parameters, according to their sensitivity indices. Temporal profiles of first order effect sensitivity indices, total sensitivity indices and those due to interactions with other model parameters have been calculated for parameters in a wastewater stabilisation pond system model. Numerical results show that the higher computational cost of global sensitivity analysis is thoroughly justified, where not only first order effects due to each parameter can be captured, but also due to interaction with other model parameters. 2. Methodology 2.1 Wastewater stabilisation pond model A differential algebraic equation system represents three stabilisation ponds in series for biological wastewater treatment. Differential equations include dynamic mass balance for the prevailing groups of bacteria: heterotrophic, autotrophic, fermenting, acetotrophic sulphate reducing and acetotrophic methanogenic bacteria; algae biomass, the different sources of organic matter, main nutrients involved in the process and methane emissions. Whereas algebraic equations describe pond interconnections, forcing functions and reaction rates of the processes. Kinetic expressions of the reaction rates and nominal value of stoichiometric coefficients were taken from Sah et al. (2011). Pond mass balances = Q , − Q , + Q , − Q , (1) , , = , , · , , , , · , + r , , − , , , ,∆ , − C , , · , m = 1 (2a) , , = r , , + , , , ,∆ , − C , , · , m = 2 (2b) Generation and consumption r , , = ν , · ρ , , (3) Qin,i, Qout,i, Qrain,i and Qevap,i stand for volumetric flow (m 3/d). Cin,i,j,k and Cout,i,j,k, are the concentrations of component j related to pond i at layer m; Vi,m is the pond volume i of the layer m (m 3) and ri,j,m is the reaction term obtained by the sum of the stoichiometric coefficients with ρk,i,m as the process reaction rates that take place within the pond i and layer m. Note that only facultative type pond has two layers. Configuration of the system was a result of the optimisation, where pond type was classified by the processes that occur within the lagoons into anaerobic, facultative and aerobic. Binary variables, ztype,i, were associated to the different types of processes, through big M constraints: Eq(4), and to different bounds of pond dimensions in order to define the pond type. Eq(5) ensured that there exists a pond of any type. ρ , − z , · M ≤ 0 k ∈ type, type= anaer, fac, aer (4) z , = 1 (5) Where anaer, fac and aer represent the sets of anaerobic, facultative and aerobic processes, respectively. The optimisation of pond stabilisation systems consists of minimising the total cost subject to environmental discharge regulations. The total cost is the sum of capital and operating costs, considering cost of land, construction, lining and maintenance of the system. Medri et al. (2007) proposed linear expressions for each type of cost: land, construction, revetment and maintenance. The objective function for minimising the total cost system is then described by Eq(6). min C = C , + C , + C , + C , (6) Where Ct is the total cost of the system (U$S) and n is the number of ponds. The model is formulated as a MINLP problem implemented in GAMS (Rosenthal, 2016) and solved by DICOPT, a solver based on the extensions of the outer-approximation algorithm for the equality relaxation strategy. More details of the model and solution strategy can be found in Ochoa et al. (2016a). 188 2.2 Global sensitivity analysis technique Sobol’s method is used to compute sensitivity indices. This method is based on variance decomposition, using Monte-Carlo simulation methods (Sobol' I., 2001; Saltelli and Tarantola, 2002). Given a function y=f(x,t), where y is a differential or algebraic state variable (e.g. heterotrophic bacteria concentration), x is a vector of k input parameters and t is the independent variable in differential equations, e.g., time; when all uncertain parameters xi vary under its probability density function, the uncertainty in y(x,t) can be quantified by its unconditional variance V(y,t). It can be decomposed as the sum of the variance of a conditional expected value and the expected value of a conditional variance, conditioning with respect to both xi and x-i. V(y,t) = V E(y|x ,t) + E V(y|x , t) (7) V(y,t) = V E(y|x ,t) + E V(y|x , t) (8) V and E correspond to variance and expected value operators, respectively. In Eq(7), V(y,t)i=V(E(y│xi,t)) computes the variance (over all possible realizations of parameter xi) of the conditional expected value of the state variable y under all parameters variation, except xi. It represents the expected reduction in the state variable variance that could be obtained if xi could be known or fixed. V(y,t)i is the first-order effect associated to parameter xi. The residual effect, E(y,t)i=E(V(y| xi,t)), is the expected value (over all realizations of parameter xi.) of the conditional variance of the state variable y when all parameters except xi change. It represents the average state variable variance if xi could be known or fixed. The same can be stated for Eq(8), by replacing xi for “all parameters except xi” (xi). Thus, the term V(y,t)i TOT=E(V(y| x-i,t)) computes the average state variable variance if all parameters except xi could be known or fixed. If equations (7) and (8) are normalised, the first-order sensitivity index, S(y,t)i and the total sensitivity index S(y,t)i TOT are defined as: S(y,t) = ( | , )( , ) = ( , )( , ) (9) S(y,t) = ( | , )( , ) = ( , )( , ) (10) To calculate sensitivity indices it is necessary to compute the unconditional and conditional variances of each state variable, involving the calculation of multiple integrals. However, Sobol’ proposes a methodology to compute the variances considering only evaluations of functions (y=f(x,t)), as defined by Eqs (11), (12) and (13) (Sobol’ and Shukhman, 2007). ∑ y − y → V(y) (11) ∑ y · y − y → V(y,t) i = 1. .k (12) ∑ y − y → V(y,t) i = 1. .k (13) A, B, Ci and Di, are matrices of dimension (N x k), N the sample size used for the Monte Carlo estimate and k the number of uncertain model parameters. Each column of both A and B matrices is a sample from the distribution function of the relative parameter. Each row is an input sample, for which a model output y can be evaluated. A is considered the ‘sample’ matrix and B the ‘re-sample’ matrix. Ci is the matrix where all parameters except xi are re-sampled, whereas Di is the matrix where only xi is re-sampled. yA, yB, yCi and yDi are vectors of N model outputs values obtained when model variables are evaluated in matrices A, B, Ci and Di, respectively. Thus, the total number of model evaluation is N(k+2). First order (S(y,t)i) and total effect (S(y,t)i TOT) indices measure the effect of the variation of the parameters on the model state variables. First sensitivity indices provide the reduction on the unconditional variance of the state variable that can be obtained if xi is fixed at its true value. On the other hand, total sensitivity indices take into account the interactions among parameters, providing information on the non-additive part of the model. The interactional sensitivity index is defined by the difference between the total sensitivity index and the first order sensitivity index. S = S − S (14) A more extended explanation of the method can be found in Ochoa et al. (2016b). 189 Uncertain Parameter Description Unit Nominal Value Standard Deviation θ Temperature coefficient 1.070 0.1338 km Mass transfer coefficient between layers in facultative pond 0.022 0.0027 KNHa Saturation/Inhibition coefficient of ammonium and ammonia nitrogen,CNH, for autotrophic bacteria, Ca gN/m 3 0.200 0.0250 inxs Fraction of nitrogen in slowly biodegradable particulate COD, Cs gN/gCOD 0.040 0.0050 KNHh Saturation/Inhibition coefficient of ammonium and ammonia nitrogen,CNH, for heterotrophic bacteria, Ch gN/m 3 0.050 0.0063 insf Fraction of nitrogen in fermentable, readily biodegradable soluble COD, Cf gN/gCOD 0.030 0.0038 ba Decay rate of autotrophic bacteria, Ca d -1 0.015 0.0019 bfb Decay rate of fermenting bacteria, Cfb d -1 0.020 0.0025 basrb Decay rate of acetotrophic sulphate reducing bacteria, Casrb d -1 0.012 0.0015 bamb Decay rate of acetotrophic methanogenic bacteria, Camb d -1 0.008 0.0010 3. Numerical results With the model presented in the previous section, a GSA analysis is performed. To generate the sample matrices N=1000 scenarios were considered. Time horizon was set to the operating time (12 days). As can be seen in Figure 1 to 4, for ammonium and ammonia nitrogen (CNH) and nitrate and nitrite nitrogen (CNO) total sensitivity indices profiles, it is evident that the temperature coefficient θ is the most influential parameter within the three ponds for almost the entire operating time. This behaviour is repeated for the 19 differential variables studied, not only due to first order effect, but also through interactional effects between the temperature coefficients with the other parameters. Therefore, focusing efforts on fixing this parameter to its true value would lead to the greatest reduction in the unconditional variance of the output. Another relevant parameter is the mass transfer coefficient between layers, km. It only contributes to the variance output in the upper and lower layer of the facultative pond. For the nitrate and nitrite nitrogen concentration, CNO,3,2, it explains up to 93% of the total variance, affecting almost the entire operating time (Figure 4b). In the case of CNH,3,2, km accounts up to the 20% of the total variance with a greater influence at the end of the time horizon, as it can be seen in Figure 2b. It also contributes to the total variance of dissolved oxygen, CO,3,2, and sulphate sulphur, CSO4,3,2, concentration of the lower layer of the facultative pond. And to a lesser extent, it influences the fermentation products, CSA,3,2, and inert soluble COD concentration, CI,3,2. (a) (b) Figure 1: Total sensitivity indices profiles for ammonium and ammonia nitrogen concentration in (a) the first aerobic pond (CNH,1,1) and (b) the second aerobic pond (CNH,2,1). 0 0.2 0.4 0.6 0.8 1 1.2 1 2 3 4 5 6 7 8 9 10 11 12 S iT O T C N H ,1 ,1 Time (d) 0 0.2 0.4 0.6 0.8 1 1.2 1 2 3 4 5 6 7 8 9 10 11 12 S iT O T C N H ,2 ,1 Time (d) θ km inxs insf KNHh KNHa ba bfb basrb bamb 190 Table 1: Uncertain parameters of the model (a) (b) Figure 2: Total sensitivity indices profile for ammonium and ammonia nitrogen concentration in the facultative pond: (a) the upper layer (CNH,3,1) and (b) lower layer(CNH,3,2). (a) (b) Figure 3: Total sensitivity indices profiles for nitrate and nitrite nitrogen concentration in (a) the first aerobic pond (CNO,1,1) and (b) the second aerobic pond (CNO,2,1). (a) (b) Figure 4: Total sensitivity indices profile for nitrate and nitrite nitrogen concentration in the facultative pond: (a) the upper layer (CNO,3,1) and (b) lower layer(CNO,3,2). 0 0.2 0.4 0.6 0.8 1 1.2 1 2 3 4 5 6 7 8 9 10 11 12 S iT O T C N H ,3 ,1 Time (d) 0 0.2 0.4 0.6 0.8 1 1.2 1 2 3 4 5 6 7 8 9 10 11 12 S iT O T C N H ,3 ,2 Time (d) θ km inxs insf KNHh KNHa ba bfb basrb bamb 0 0.2 0.4 0.6 0.8 1 1.2 1 2 3 4 5 6 7 8 9 10 11 12 S iT O T C N O ,1 ,1 Time (d) 0 0.2 0.4 0.6 0.8 1 1.2 1 2 3 4 5 6 7 8 9 10 11 12 S iT O T C N O ,2 ,1 Time (d) θ km inxs insf KNHh KNHa ba bfb basrb bamb 0 0.2 0.4 0.6 0.8 1 1.2 1 2 3 4 5 6 7 8 9 10 11 12 S iT O T C N O ,3 ,1 Time (d) 0 0.2 0.4 0.6 0.8 1 1.2 1 2 3 4 5 6 7 8 9 10 11 12 S iT O T C N O ,3 ,2 Time (d) θ km inxs insf KNHh KNHa ba bfb basrb bamb 191 Saturation/Inhibition coefficient of CNH for Ca, KNHha, accounts for 65% and 95% of the total variance of CNH at the beginning of the process in the second aerobic pond and in the upper layer of the facultative pond, respectively (Figure 1b and Figure 2a). This fact is in concordance with the high growth rate of the autotrophic bacteria at that point of the time horizon (Ca,2,1 and Ca,3,1). This parameter also contributes to the total variance of nitrate and nitrite nitrogen of the second aerobic pond, CNO,2,1, (Figure 3b) through its first order and interactional effect, and to the total variance of nitrate and dissolved oxygen in the second aerobic pond and in the upper layer of the last pond (CO,2,1 and CO,3,1) only by its interactional effects. The influence is again appreciated at the beginning of the operating time. Finally, the variance output is comprised of two more parameter: the Saturation/Inhibition coefficient of CNH for Ch, KNHh, and the fraction of nitrogen in Cs, inxs, and. The contribution of both parameters is low; the former affects the fermentation products (CSA,2,1 and CSA,3,1), fermentable, readily biodegradable soluble COD (CF,2,1 and CF,3,1) and dissolved oxygen (CO,1,1) concentration profiles. Whereas, the latter only influences the nitrate and nitrite nitrogen variance output less than 18% along all operating days, as it can be appreciated in Figure 3b and Figure 4a. 4. Conclusions A proper knowledge of the influence of the model parameters on the state variables has allowed their classification and provides useful information for parameter estimation. GSA allows determining which parameters are the most influential through their effects due to interactions with other parameters and due to first order effects along the entire time horizon. The temperature coefficient θ is the parameter that has the strongest first order influence during the time horizon in the wastewater stabilisation pond system model, which indicates the need of its proper determination. It is important to estimate the most accurate values of other parameters like km and KNHha in order to avoid their first order effects that explain 100% of the variance during some intervals within the time horizon of some variables output. Such information is crucial for identifying the set of these three parameters that need to be determined more precisely, as stated in the results section, and allows the identification of the set of parameters which has no contribution to the total output variance. Acknowledgments Authors greatly acknowledge the National Research Council CONICET, ANPCyT and Universidad Nacional del Sur for supporting their work through grants PIP 2011 11220110101078 (2013-2015), PICT 2012-2469 and PGI 24/M125 respectively. References Medri, W., Costa, R.H.R., Medri, V., Belli Filho, P., 2007, Stabilization Pond Systems: Cost Estimation for the Treatment of Piggery Waste, Transactions of the American Society of Agricultural and Biological Engineers, 50(4):1409-1414. Ochoa, M.P., Estrada, V., Hoch, P.M., 2016a, MINLP Wastewater Stabilisation Ponds Synthesis using Rigorous Models under Different Scenarios, Computer Aided Chemical Engineering. Ochoa, M.P., Estrada, V., Di Maggio, J., Hoch, P.M., 2016b, Dynamic global sensitivity analysis in bioreactor networks for bioethanol production, Bioresource Technology, 200: 666-679. PSEnterprise Ltd., 2014, gPROMS Model Developer Guide. London: Process Systems Enterprise Ltd. Rosenthal, R.E., 2016, GAMS- A User’s Guide, Washington, DC, USA. Sah, L., Rousseau, D. P. L., Hooijmans, C.M., Lens, P.N.L., 2011, 3D model for a secondary facultative pond. Ecological Modelling, 222(9), 1592-1603. Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., Tarantola, S., 2008, Global Sensitivity Analysis. The Primer. John Wiley & Sons, Ltd. Saltelli, A., Tarantola, S. (2002). On the relative importance of input factors in mathematical models: safety assessment for nuclear waste disposal. Journal of the American Statistical Association, 97, 702–709. Sobol', I. M. (1993). Sensitivity estimates for nonlinear mathematical models. Mathematical Modelling & Computational Experiment, 1, 407-414. Sobol', I. (2001). Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates. Mathematics & Computers in Simulation, 55, 271-280. Sobol’, I., Shukhman, B. (2007). On global sensitivity indices: Monte Carlo estimates affected by random errors. Monte Carlo Methods and Appl., 13(1), 89-97. 192