Accounting for erroneous model structures in biokinetic process models

Accounting for erroneous model structures in biokinetic process models Kris Villez a,b , Dario Del Giudice a,c,d , Marc B. Neumann e,f , Jörg Rieckermann a a Eawag: Swiss Federal Institute of Aquatic Science and Technology, Überlandstrasse 133, 8600 Dübendorf, Switzerland b ORNL: Oak Ridge National Laboratory, Oak Ridge, TN, USA c ETHZ: Swiss Federal Institute of Technology, Stefano-Franscini-Platz 5, 8093 Zürich, Switzerland d Department of Civil, Construction & Environmental Engineering, NC State University, Mann Hall 311, 2501 Stinson Drive, Raleigh, NC, 27695, USA e Basque Centre for Climate Change (BC3), Scientic Campus of the University of the Basque Country, Sede Building 1, 1st oor, 48940 Leioa, Spain f IKERBASQUE, Basque Foundation for Science, Maria Diaz de Haro 3, 6 solairua, 48013 Bilbao, Spain Abstract In engineering practice, model-based design requires not only a good processbased model, but also a good description of stochastic disturbances and measurement errors to learn credible parameter values from observations. However, typical methods use Gaussian error models, which often cannot describe the complex temporal patterns of residuals. Consequently, this results in overcondence in the identied parameters and, in turn, optimistic reactor designs. In this work, we assess the strengths and weaknesses of a method to statistically describe these patterns with autocorrelated error models. This method produces increased widths of the credible prediction intervals following the inclusion of the bias term, in turn leading to more conservative design choices. However, we also show that the augmented error model is not a universal tool, as its application cannot guarantee the desired reliability of the resulting wastewater reactor design. Keywords: bias description; kinetic model; process design; wastewater treatment; uncertainty Email address: [email protected] (Kris Villez) Preprint submitted to Reliability Engineering and System Safety July 21, 2020 This document is the Accepted Manuscript version of a Published Work that appeared in final form in: Villez K., Del Giudice D., Neumann M.B., Rieckermann J. 2020. Accounting for erroneous model structures in biokinetic process models. RELIABILITY ENGINEERING & SYSTEM SAFETY. 203. DOI (10.1016/ j.ress.2020.107075). © 2020 Elsevier Ltd This manuscript version is made available under the CC-BY-NC-ND 3.0 license http://creativecommons.org/licenses/ by-nc-nd/3.0/ Copyright notice This manuscript has been authored in part by UT-Battelle, LLC, under contract DE-AC05-00OR22725 with the US Department of Energy (DOE). The US government retains and the publisher, by accepting the article for publication, acknowledges that the US government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for US government purposes. DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doepublicaccess-plan). 1. Introduction 1 In current environmental engineering practice, deterministic process-based 2 modeling is a common tool to better understand the functioning of complex 3 wastewater collection and treatment systems. The gold standard is to im4 prove prediction performance of our models by tting them to observations. 5 Consequently, the advent of ubiquitous sensing leads to an unintended yet 6 commonly observed situation where sensors reveal more details than mech7 anistic models can capture. When this is the case, uncertainty estimates 8 obtained from statistical inference with mechanistic models are almost cer9 tainly too narrow as the applied model structure is too restrictive relative 10 to the observed reality. A long-standing question is therefore whether risk11 based design, based on uncertainty estimates from statistical inference with 12 mechanistic models, is actually feasible. In this work, we test one method 13 designed to address this issue to a case of WWTP design and discuss its 14 potential and limitations. 15 Accounting for model parameter uncertainty is crucial for risk-based decision16 making, including infrastructure design and operations (e.g., Cagno et al., 17 2011;Scheidegger et al.,2013;Kabir et al.,2015;Scheidegger et al.,2015; 18 Jensen and Jerez,2018). Conventional methods for uncertainty analysis are 19 based on a two-step approach, consisting of (a) quantication of input un20 certainty, measurement uncertainty, and subsequent uncertainty of model 21 parameters followed by (b) propagation of the quantied uncertainty to the 22 system performance measure of interest (e.g., Van Griensven and Meixner, 23 2007;Sin et al.,2009;Guo and Murphy,2012;Del Giudice et al.,2016). 24 2 However, it has been demonstrated before how systematic deciencies in 25 model structure, next to input and measurement uncertainty, also lead to bi26 ased model parameters and, consequently, incorrect design of infrastructural 27 elements, such as biological reactor systems (Neumann and Gujer,2008). 28 Unfortunately, to our knowledge, no one has attempted to provide a method 29 to solve this particular problem, i.e. to identify systematic discrepancies in 30 process-based models so to account for them during model-based design. 31 Recently, statisticians have been suggesting a promising approach to solve 32 this dilemma. The underlying idea is to not assume identically and indepen33 dently distributed (i.i.d.) errors for mismatches between models and obser34 vations (Liu and Zachara,2001), but to explicitly account for mismatches by 35 adding a stochastic auto-correlated process to the i.i.d. measurement error 36 model (Craig et al.,2001;Kennedy and O'Hagan,2001;Bayarri et al.,2007). 37 This is known as the bias description method and focuses on the modeling 38 of the symptoms of a mismatch between model structure and reality. While 39 this does not identify or tackle the root cause of these symptoms, it has 40 been proven to be a computationally ecient tool to increase the reliability 41 of model-based predictions compared to standard regression approaches in a 42 variety of systems from lakes to natural catchments to urban hydrology (Di43 etzel and Reichert,2012;Reichert and Schuwirth,2012;Del Giudice et al., 44 2015). Therefore, we expect that the bias description method also improves 45 the reliability of predictions with structurally decient wastewater treatment 46 models in view of risk-based design. Specically, adding a stochastic mea47 surement error term to a model given the same amount of experimental 48 measurements is expected to reduce the relative information-richness of the 49 experimental data and lead to larger credibility intervals of model parame50 ters, wider prediction intervals and, by avoiding overcondent predictions, a 51 more trustworthy design. 52 Note that the bias description method can be regarded as a grey box or 53 hybrid modelling strategy. Indeed, the resulting model consists of a mech54 anistic model for the studied process (white box) and a stochastic model 55 for auto-correlated measurement errors (black box). Other grey box ap56 proaches may be based on the inclusion of time-variant parameters (Reichert 57 and Mieleitner,2009;Lin and Beck,2012) or integration of non-parametric 58 elements into a model structure that is mechanistic otherwise (Ma²i¢ et al., 59 2017). 60 In this contribution, we apply the bias description method to investigate 61 the impact of model structure decits for process design. We use Neumann 62 3 and Gujer (2008) as a benchmark to evaluate the benets and limitations 63 of the bias description method for model-based design and refer to it as the 64 reference study . While this reference study concerns a conceptually simple 65 case, using it in this study highlights (a) that the apparent simplicity of this 66 case is rather deceptive and (b) that challenges associated with model-reality 67 mismatch are to be expected for both simple and complex systems. 68 2. Material and methods 69 2.1. Applied error models 70 In a vast majority of environmental modeling studies, the measurement 71 error is assumed to be i.i.d. For example, Hauduc et al. (2015) compares an 72 extensive list of model performance criteria for wastewater treatment mod73 elling yet does not list any criterion which accounts for autocorrelated model 74 prediction errors. One approach considered in Cierkens et al. (2012), consists 75 of downsampling time series to avoid the appearance of autocorrelation. As 76 explained in the same study, this leads to an inecient use of the available 77 data and, more importantly, cannot account at all for model structure decits 78 as a potential root cause of autocorrelated residuals. Ignoring the presence 79 of autocorrelated residuals was shown to lead to overcondence in the pro80 duced model and, subsequently, poor decision-making, as was shown also in 81 the reference study. Most often, a Gaussian distribution is assumed for the 82 measurement errors. Such a model of measurement error cannot account for 83 systematic deviations between the assumed model and the observed mea84 surements, i.e. bias. One way of accounting for bias is by adding terms, such 85 as a stochastic autocorrelated error bias term, to the measurement equation 86 (Craig et al.,2001;Kennedy and O'Hagan,2001;Bayarri et al.,2007). In 87 this work, we describe the observable output time-series (i.e., measured con88 centration, yo ) as a sum of a deterministic dynamic model output ( y , the 89 modeled concentration), a classical Gaussian measurement error ( e(ψ) ), and 90 an auto-correlated error term ( b(ψ) ): 91 yo(θ, γ, ψ) = y(θ) + γ+ b(ψ) + e(ψ) (1) where θ and γ are parameters of the deterministic parts of the model and 92 ψ are those of the stochastic parts (errors). The bias term b(ψ) decribes an 93 autocorrelated error ( b(ψ)∼ N(0,Σb(σb, τ)) ) and can be included to account 94 4 for time-dependent deviations between model and observations (Reichert and 95 Schuwirth,2012). Note that this bias term represents a stochastic process, 96 thus describing aleatory uncertainty, although the deviations between model 97 and observations may actually be systematic, possibly even deterministic. 98 These deviations are expected to be systematic when they are caused by a 99 lack of knowledge about the true data-generating process. This lack of knowl100 edge is typically characterized as a source of epistemic uncertainty rather than 101 aleatory uncertainty. 102 The bias term has two parameters, the standard deviation σb and the 103 correlation length τ : 104 Σb(i, j) := σb2·e−|ti−tj|2/τ . (2) The random measurement error is temporally independent ( e∼ N(0,Σe(σe)) ) 105 and is characterized by the parameter σe : 106 Σe(i, j) := (σe2, i =j 0, i 6=j. (3) Together, these error terms with parameters ψ={τ, σb, σe} account for 107 the fact that the deterministic model may not reproduce the modeled data 108 set exactly. Note that the symbols σb and σe are chosen to convey the idea 109 that they both describe the magnitude of variation of a stochastic term in the 110 measurement equation. The symbol for the correlation length, τ , is chosen 111 to highlight the fact that it describes a time-scale. 112 The statistical formulation in Eq. 1naturally leads to the likelihood 113 function L(yo|θ, ψ) which describes how likely the considered model with 114 parameters (θ, ψ) generated the recorded data, yo . The likelihood of the 115 measurements conditional to the model parameters is: 116 L(yo|θ, ψ) = (2π)−n 2 qdetΣ exp −1 2hyo−yiT (Σ)−1hyo−yi (4) where Σ is the variance-covariance matrix for the stochastic deviations 117 between model and observations: 118 5 Σ := Σe+ Σb, (5) with Σb and Σe dened in 2and 3, and with n equal to the number of 119 measurements. In order to interpret the results of parameter estimation, we 120 dene the parameters α and σ such that σ2 e:= (1 −α)σ2 and σ2 b:= α σ2 . 121 This means we can express the variance-covariance matrix above equivalently 122 as: 123 Σ := σ2·h(1 −α)In×n+α Ki (6) K(i, j) := e−|ti−tj|2/τ (7) In this form, σ is a measure for the overall spread of the deviations be124 tween the model and the measurements and α is a parameter that denes 125 the relative importance of the bias in the overall variance-covariance matrix. 126 Meaningful values for α are between 0 and 1, with α= 1 leading to the omis127 sion of the independent measurement noise ( σe= 0 ) and α= 0 expressing 128 that there is no bias ( σb= 0 ). Note that setting α= 0 reproduces the model 129 without a bias term. Put otherwise, the model with bias term includes the 130 model without bias term as special case. 131 The addition of a bias term accounts for underestimation of parameter 132 uncertainty when a conventional yet unrealistic distribution for the model 133 error is assumed (e.g., uncorrelated). This is expected to produce a wider 134 predictive distribution, possibly leading to a better quantication of and a 135 reduction of the risk of under-design or over-design. However, special atten136 tion must be given to the parameter estimation method as increased model 137 exibility can lead to unidentiability (see e.g., Renard et al.,2010). 138 2.1.1. Model parameter estimation 139 We apply a Bayesian approach for two reasons. First, we favour a Bayesian 140 framework as a way to make prior beliefs explicit. Second, without any form 141 of prior, some of the parameters of the variance-covariance matrix Σ can be 142 structurally unidentiable (for denitions, see Dochain et al.,1995;Dochain 143 and Vanrolleghem,2001;Petersen et al.,2003). More specically, when τ= 0 144 the matrix K equals the identity matrix and likelihood L(yo|θ, ψ) becomes 145 insensitive to the value of α . As a result, no unique value for α can be 146 6 identied under any circumstances as long as τ= 0 , i.e. α is structurally 147 unidentiable. In the formulation with σe and σb , any increase of σe can 148 be compensated exactly by an equivalent decrease of σb when τ= 0 . For 149 small values of τ , e.g. close to the measurement interval or smaller, this is 150 expected to lead to a lack of practical identiability, even if structural identi151 ability could be guaranteed in principle. In early experiments with uniform 152 priors for τ , we observed that this can induce a lack of convergence and poor 153 mixing conditions for the applied sampling methods, similar to observations 154 described in Renard et al. (2010). Applying an informative prior solves this 155 identiability problem and can therefore also be interpreted as a form of 156 regularization (e.g., Scales and Tenorio,2001;Murphy,2012;Hastie et al., 157 2015). 158 Bayesian calibration aims at characterizing the distribution described by 159 the posterior likelihood L(θ, ψ|yo)∝ L(yo|θ, ψ)· L(θ, ψ) , where the prior 160 likelihood L(θ, ψ) expresses the prior beliefs about the parameters. In this 161 work, the posterior distribution is approximated with a sample of L(θ, ψ|yo) 162 drawn with a Markov Chain Monte Carlo (MCMC) sampler (see Numerical 163 implementation). 164 2.2. Biokinetic model parameter identication with batch experiments 165 To study the eects of model structure error and the utility of the bias de166 scription method, we execute simulations with the dynamic biokinetic model 167 used in the reference study. Concretely, a series of batch experiments is sim168 ulated in which a substrate, with concentration s(t) , is consumed by a cell 169 culture with a xed concentration. The conversion rate r(t) depends on the 170 substrate by means of time-invariant Tessier kinetics so that one can write: 171 ds(t) dt =−r(t) (8) s(t= 0) = s0 (9) r(t) = rTessier max ·1−exp −s(t) KTessier S·x(t) (10) During the experiment, noisy measurements of the true substrate con172 centration, yo(t) , are simulated by the following measurement error model, 173 which is a zero-mean Gaussian noise term: 174 7 yo(t) = s(t) + e(t) (11) e(t)∼N(0, σe) (12) Fixed parameters for each simulation are the same as in the reference 175 study: s0 (initial substrate concentration, 5g/m3 ), rmax (maximum con176 version rate, 1g/m3.h ). The simulated time is T= 8 hours. The anity 177 constant ( KTessier S ) and measurement error standard deviation ( σe ) are varied 178 yet constant in every simulated experiment. KTessier S is varied from 0.1g/m3 179 to 1.5g/m3 in steps of 0.2g/m3 . This allows simulating a wide range of 180 process conditions, including both low and high values for KTessier S relative 181 to the initial substrate concentration. Two values for the simulated σe are 182 considered, as in the reference study. In the low-noise case, σe takes the value 183 0.01 g/m3 . In the high-noise case, it takes the value 0.1g/m3 . The vector 184 θ equals s0, rmax, KT essier S T . The two simulated noisy time series obtained 185 with KTessier S= 0.7g/m3 are shown in the supplementary information (Fig. 186 S.1 and Fig. S.2). 187 For each simulation experiment, parameter identication is executed with 188 four distinct model structures. The rst model matches the above model 189 structure (Eq. 8-Eq. 12) exactly. This represents an idealized situation where 190 the structure of the calibrated model matches reality (ground truth) exactly. 191 The identied parameters are S0 , µTessier max , KTessier S , and σe . A second model 192 is obtained by replacing the Tessier kinetics with the alternative and more 193 commonly used Monod kinetics. Practically, Eq. 10, is replaced with the 194 following equation: 195 r(t) = rMonod max ·s(t) KMonod s+s(t) (13) The estimated parameters are now s0 , µMonod max , KMonod s , and σe with 196 θ=s0, rmax, KMonod S T . This case represents the likely situation that a 197 modeling practitioner uses the common-place Monod model structure and 198 does not observe the bias that results. This is very likely in the high-noise 199 case (see reference study). Given this diculty, the stochastic bias term de200 scribed above is included to capture the systematic deviations between the 201 model predictions and measurements. To achieve this, the previously applied 202 measurement equation (Eq. 11) is replaced with the following equations: 203 8 yo(t) = s(t) + b(t) + e(t) (14) e(t)∼ N(0, σe) (15) b∼ N(0,Σb(σb, τ)) (16) with the parameters τ dened as above and σb and σe reparametrized with 204 α and σ . This results in a third model, where a Tessier model is combined 205 with the statistical bias description and which requires specication of the 206 parameters s0 , µTessier max , KTessier S , σ , α , and τ . As the Tessier model has 207 the same structure as the data-generating model, one can expect a good 208 model t with α close to 0 and estimates of s0 , µTessier max , KTessier S , and σ 209 that are close to ground truth values. Finally, the fourth model combines 210 the presumed Monod kinetics with the statistical bias description and the 211 identied parameters are s0 , µMonod max , KMonod s , σ , α , and τ . In this case we 212 can expect that the present model structure bias is accommodated by means 213 of the statistical bias description. If so, this should increase the width of the 214 prediction intervals and thereby improve the reliability of the model (Reichert 215 and Schuwirth,2012). 216 2.3. Numerical implementation 217 The biokinetic model, parameter estimation, and uncertainty propaga218 tion were implemented in Matlab (R2019a). The prior probabilities for the 219 parameters were set based on the authors' experience. They are all indepen220 dent of each other. All priors are uniform, except for σ and τ . The prior 221 likelihood for σ is proportional to its inverse and is equivalent to the Jereys 222 prior conditional to xed values for all other parameters (see Box and Tiao, 223 1973). The prior likelihood for τ is the sine function supported between 0 and 224 2T . This prior equals zero at τ= 0 and τ= 2 T and one at τ=T . This ex225 presses the subjective belief that the autocorrelation length of the deviations 226 due to model structure error is expected to be similar to the duration of the 227 experiment. The priors are specied completely in Table 1. We rst run an 228 adaptive MCMC algorithm (Vihola,2012) to nd a good guess for the maxi229 mum a posteriori estimates and a good proposal variance-covariance matrix. 230 With these results, we execute a (non-adaptive) MCMC algorithm to obtain 231 20,000 samples from L(θ, ψ|yo) . The rst 10,000 samples are considered to 232 correspond to the burn-in phase of the sampler, during which eects of the 233 initial sample may still be apparent. These samples are therefore discarded, 234 as is common in practice (Gilks et al.,1996). 235 9 Figure 4: Distributions of the ratio of predicted steady state concentrations to the ground truth concentration as a function of the anity constant ( KS )  High noise case ( σe= 0.1g/m3 ). Red horizontal whiskers indicate the two-sided 99% credible intervals. Left side beans: without bias description; Right side beans: with bias description. Top: Tessier model - All credible intervals include the ideal ratio (equal to 1), except for the simulation with KS= 0.1g/m3 without bias term. The uncertainty increases when a bias term is added to the model. Bottom: Monod model  Including the bias term in the model increases the reliability of the credible intervals. These intervals include the ideal ratio for two cases ( KS= 0.5 and 0.7g/m3 ). For KS= 0.3g/m3 this already amounts to 8.6%. Since the model structure 354 error primarily relates to the curvature of the conversion rate in this region, 355 it follows that model structure error will always be dicult to detect when 356 this time fraction is low. 357 In the supplementary information, we provide results obtained with the 358 modied method. We omit the information obtained during the washout ex359 periment during prediction. In this case, the uncertainty in the predictions 360 16 Figure 5: Distributions of the parameter α for both models with a bias term in all highnoise cases ( σe= 0.1g/m3 ). When the Tessier model is selected (no model structure error), the posterior distribution of α is shifted to the left of the prior, thus suggesting the kinetic model structure is adequate. In contrast, the posterior of α is similar to or located at the right of the prior when the Monod model is used in all but one case ( KS= 0.1g/m3 ), thus providing a useful indication of model structure error. is reduced signicantly to the point that none of the 99% credible intervals 361 include the ground truth (see Fig.S.3). This is explained by the fact that 362 the estimates for rmax and KS exhibit strong correlation (see supplemen363 tary information for details). However, since the modication relates to the 364 prediction step only, one can still use the posteriors for α as a detection 365 mechanism for bias. 366 17 4. Discussion 367 4.1. Summary and limitations of the experimental simulation study 368 Summary. The numerical results described above suggest that the inclusion 369 of an additive auto-correlated error process into a measurement error model 370 can improve the reliability of model-based designs. This is true even when 371 only a subset of the identied parameters are used during prediction (here we 372 only used the estimates for KS ) and even when the experimental setting for 373 prediction (steady state) is dierent from the experimental conditions used 374 for model identication (batch experiment). In our case, the bias description 375 method improves the reliability in all cases. Despite this improvement, the 376 computed credible intervals include the ground truth value only in a lim377 ited number of cases with model structure error, meaning that guaranteed 378 reliability cannot be obtained with the studied method. Thus, the inclusion 379 of a bias description term for the purpose of prediction can be advised as 380 a relatively fast and easy way to account for errors in the proposed model 381 structure, however only when one is unable to modify the model structure 382 itself. This is especially relevant in engineering applications where one is re383 stricted to specic process representations (e.g., Monod kinetics) or software 384 with limited exibility. While the bias description method improves the re385 liability of the model predicitions only in a limited way, it is very useful as 386 a tool to detect the presence of bias during model identication, especially 387 when reformulated with the α parameter. 388 Limitations. In this study, a simple case was chosen deliberately for two 389 reasons. First, this enabled an objective comparison of the bias descrip390 tive method with the historical results in the reference study (Neumann and 391 Gujer,2008). Second, the apparent simplicity of the case also highlights the 392 challenge of generating reliable predictions with mechanistic models, induced 393 by the typical lack of exibility of such models. The chosen scope also means 394 that our study comes with some limitations, which are: 395  The general applicability of the bias description method is not demon396 strated. However, the bias description method could easily be adapted 397 to more complex systems. One could incorporate a bias term to ex398 press correlation between multiple measurements, of the same or dis399 tinct variables measured in the same location or dierent locations. In 400 this case, the covariance between two measurements, as expressed by Σ , 401 would not only be a function of (a) the time dierence ( ti−tj , see (6)), 402 18 as in our study, but also of (b) spatial distance in one or more dimen403 sions and (c) eects of measurement error correlation between distinct 404 sensors measuring the same or distinct variables. This generalization 405 of the present model is likely most convenient when the bias error term 406 is modelled as spatio-temporal Gaussian process (e.g., De Cesare et al., 407 2001;Gneiting,2002;Stein,2005). 408  The methods applied in both the reference study and ours are based on 409 methods that account for aleatory uncertainty only. However, the lack 410 of knowledge about the model structure is typically epistemic in nature 411 and may therefore be dicult to account for in this way. Epistemic 412 uncertainty may however be reduced by using more exible models 413 (Ma²i¢ et al.,2017) while increasing parametric uncertainty, which can 414 be handled as an aleatory source of uncertainty with currently available 415 methods. Still, the adoption of alternative frameworks for uncertainty 416 analysis (Parsons,2001;Rao et al.,2008) may be suited to handle 417 epistemic uncertainty directly. In summary, the handling of epistemic 418 uncertainty deserves more attention. 419 4.2. General consequences for practical uncertainty and reliability analysis 420 Utility of the bias description method. In our opinion, the detection of sys421 tematic deviations between the assumed model structure and the data-generating 422 process is the most useful feature of the bias description method. For this 423 reason, we recommend that a model is inspected for bias by (a) adding a bias 424 term in the assumed model, specifying a prior for alpha concentrated around 425 a strictly positive value, as suggested here, and (b) inspecting the posterior 426 of α whenever an inappropriate model structure is suspected. Reformulation 427 of the error model (bias + measurement error) with α , σ , and τ proved very 428 helpful as it enables interpreting α as an indicator for the relative impor429 tance of model structure error. In cases where the posterior probability mass 430 is not shifted towards zero, relative to the prior, the modeler should suspect 431 the presence of bias. When this is detected, potential model improvements 432 may include the use of time-dependent parameters (Reichert and Mieleitner, 433 2009;Lin and Beck,2012) and/or input errors (Del Giudice et al.,2016) or 434 a change in model structure (Del Giudice et al.,2015;Ma²i¢ et al.,2017). 435 While the method increases the reliability of the obtained steady-state pol436 lutant concentration predictions, it is important to note that the observation 437 of this benet depends strongly on the root cause of the observed bias. For 438 19 this reason, detection of bias should be followed by exploratory analysis of 439 the residuals and development of a better model structure (e.g., Reichert 440 and Mieleitner,2009;Del Giudice et al.,2013). We do not recommend ex441 ploiting the bias term for prediction without search for the underlying causes 442 for model structure decits, especially considering that the ground truth is 443 rarely included in the produced credible intervals. Ultimately, the utility 444 of any approach depends on whether it can successfully describe the rele445 vant sources of the deviations between model predictions and the measured 446 variables (Brynjarsdóttir and O'Hagan,2014;Wani et al.,2019). 447 Parameter interpretation and transferability. The mechanistic interpretation 448 of identied values for the parameters in the deterministic part of the model 449 is nearly impossible when bias is present. Adding an auto-correlated additive 450 error term contributes to a better reliability of the model predictions but can451 not provide a clearer interpretation of the parameter values or a direction to 452 a more appropriate model structure. Indeed, the parameter estimates remain 453 biased. Importantly, this is a likely scenario in wastewater engineering due to 454 the extremely simplied representation of biological processes during model 455 construction. Furthermore, obtaining proofs of a lack of bias is extremely 456 dicult to achieve so that a straightforward interpretation of parameter val457 ues is unlikely, even when the model structure may be appropriate. However, 458 grey-box or hybrid models may oer intepretability and transparency at the 459 cost of computational eorts (see introduction above). 460 Data quality. The quality of the simulated measurements in the studied case 461 is fairly high relative to current experience in the wastewater sector. However, 462 sensor hardware has become increasingly robust in the last three decades 463 (Olsson,2012) and there is no obvious reason why this trend should stop 464 now. It is therefore reasonable to expect that the presence of bias can be 465 detected easily in the future, either by statistical tests for auto-correlation 466 of the residuals, as in the reference study, or with descriptive methods, as in 467 this study. This will also facilitate the modication of the model structure 468 in accordance to the envisioned high-quality data. 469 4.3. Future work 470 Through this work, we identied several avenues of further research. 471 These include: 472 20  Develop and study methods for parameter estimation and parameter 473 interpretation under presence of model structure error. 474  Develop a systematic approach to the formulation of prior distributions, 475 especially when exibility is at odds with model structure or parameter 476 identiability. 477  Evaluation of experimental design methods to improve the chances of 478 detection of model structure errors. 479  Adopt and evaluate methods to handle epistemic uncertainty in model480 based process design and operation. 481 5. Conclusions 482 In this paper, we investigated the challenge of structural model decits in 483 risk-based reactor design. This is a relevant problem, because digitalization 484 will improve sensor resolution and spatial coverage of reactors, which will 485 reveal mismatches (i.e, bias) in our common engineering models (which have 486 been developed in the data-scarce past, often by grab sampling). Auto487 correlated mathematical formulations have been suggested to improve the 488 description of such biases. 489 In summary, our study shows that 490  Adding auto-correlation terms in the measurement error model as a 491 way to account for model structure decits signicantly improves the 492 reliability of biokinetic models. 493  Bias description enables accounting for predictive uncertainty during 494 process design to a large degree. This does not produce a guaranteed 495 reliability of the resulting design however. It is therefore not a bullet496 proof solution to the presence of model-reality mismatch. 497  The studied bias description method is an adequate tool to identify the 498 presence of model structure decits in presence of noisy experimental 499 data. 500 21 Acknowledgments 501 We thank Peter Reichert and Sanda Dejanic for their helpful insight into 502 the studied problem. Marc B. Neumann acknowledges nancial support pro503 vided by the Spanish Government through the BC3 María de Maeztu ex504 cellence accreditation 2018-2022 (MDM-2017-0714) and the Ramón y Cajal 505 grant (RYC-2013-13628); and by the Basque Government through the BERC 506 2018-2021 program. 507 6. Author contributions 508 Dario del Giudice: Methodology, Software, Data analysis, Writing - Re509 view & Editing. Marc B. Neumann: Methodology, Software, Data analysis, 510 Writing - Review & Editing. Jörg Rieckermann: Conceptualization, Method511 ology, Writing - Review & Editing. Kris Villez: Methodology, Software, 512 Formal analysis, Data analysis, Visualization, Writing - Original Draft. 513 References 514 Bayarri, M.J., Berger, J.O., Paulo, R., Sacks, J., Cafeo, J.A., Cavendish, J., 515 Lin, C.H., Tu, J., 2007. A framework for validation of computer models. 516 Technometrics 49, 138154. doi: 10.1198/004017007000000092 . 517 Box, G.E.P., Tiao, G.C., 1973. Bayesian Inference in Statistical Analysis. 518 Addison-Wesley, Reading, MA, USA. 519 Brynjarsdóttir, J., O'Hagan, A., 2014. Learning about physical parame520 ters: The importance of model discrepancy. Inverse Problems 30, 114007. 521 doi: 10.1088/0266-5611/30/11/114007 . 522 Cagno, E., De Ambroggi, M., Grande, O., Trucco, P., 2011. Risk analysis 523 of underground infrastructures in urban areas. Reliability Engineering & 524 System Safety 96, 139148. doi: 10.1016/j.ress.2010.07.011 . 525 Cierkens, K., Plano, S., Benedetti, L., Weijers, S., de Jonge, J., Nopens, I., 526 2012. Impact of inuent data frequency and model structure on the quality 527 of WWTP model calibration and uncertainty. Water Science & Technology 528 65, 233242. doi: 10.2166/wst.2012.081 . 529 22 Craig, P.S., Goldstein, M., Rougier, J.C., Seheult, A.H., 2001. Bayesian 530 forecasting for complex systems using computer simulators. Journal 531 of the American Statistical Association 96, 717729. doi: 10.1198/ 532 016214501753168370 . 533 De Cesare, L., Myers, D.E., Posa, D., 2001. Estimating and modeling space534 time correlation structures. Statistics & Probability Letters 51, 914. 535 doi: 10.1016/S0167-7152(00)00131-0 . 536 Del Giudice, D., Albert, C., Rieckermann, J., Reichert, P., 2016. Describing 537 the catchment-averaged precipitation as a stochastic process improves pa538 rameter and input estimation. Water Resources Research 52, 31623186. 539 doi: 10.1002/2015WR017871 . 540 Del Giudice, D., Honti, M., Scheidegger, A., Albert, C., Reichert, P., Rieck541 ermann, J., 2013. Improving uncertainty estimation in urban hydrological 542 modeling by statistically describing bias. Hydrology and Earth System 543 Sciences 17, 42094225. doi: 10.5194/hess-17-4209-2013 . 544 Del Giudice, D., Reichert, P., Bares, V., Albert, C., Rieckermann, J., 2015. 545 Model bias and complexity - Understanding the eects of structural decits 546 and input errors on runo predictions. Environmental Modelling & Soft547 ware 64, 205214. doi: 10.1016/j.envsoft.2014.11.006 . 548 Dietzel, A., Reichert, P., 2012. Calibration of computationally demand549 ing and structurally uncertain models with an application to a lake wa550 ter quality model. Environmental Modelling & Software 38, 129146. 551 doi: http://dx.doi.org/10.1016/j.envsoft.2012.05.007 . 552 Dochain, D., Vanrolleghem, P.A., 2001. Dynamical modelling & estimation 553 in wastewater treatment processes. IWA Publishing, London, UK. 554 Dochain, D., Vanrolleghem, P.A., Van Daele, M., 1995. Structural identia555 bility of biokinetic models of activated sludge respiration. Water Research 556 29, 25712578. doi: 10.1016/0043-1354(95)00106-U . 557 Gilks, W.R., Richardson, S., Spiegelhalter, D.J., 1996. Markov Chain Monte 558 Carlo in practice. Chapman and Hall, London, UK. 559 23 Gneiting, T., 2002. Nonseparable, stationary covariance functions for space560 time data. Journal of the American Statistical Association 97, 590600. 561 doi: 10.1198/016214502760047113 . 562 Guo, M., Murphy, R.J., 2012. LCA data quality: sensitivity and uncertainty 563 analysis. Science of the Total Environment 435, 230243. doi: 10.1016/j. 564 scitotenv.2012.07.006 . 565 Hastie, T., Tibshirani, R., Wainwright, M., 2015. Statistical learning with 566 sparsity: the Lasso and generalizations. Chapman and Hall/CRC, Boca 567 Raton, FL, USA. 568 Hauduc, H., Neumann, M.B., Muschalla, D., Gamerith, V., Gillot, S., Van569 rolleghem, P.A., 2015. Eciency criteria for environmental model quality 570 assessment: A review and its application to wastewater treatment. Envi571 ronmental Modelling & Software 68, 196204. doi: 10.1016/j.envsoft. 572 2015.02.004 . 573 Jensen, H.A., Jerez, D.J., 2018. A stochastic framework for reliability and 574 sensitivity analysis of large scale water distribution networks. Reliability 575 Engineering & System Safety 176, 8092. doi: 10.1016/j.ress.2018.04. 576 001 . 577 Kabir, G., Tesfamariam, S., Sadiq, R., 2015. Predicting water main failures 578 using Bayesian model averaging and survival modelling approach. Relia579 bility Engineering & System Safety 142, 498514. doi: 10.1016/j.ress. 580 2015.06.011 . 581 Kennedy, M.C., O'Hagan, A., 2001. Bayesian calibration of computer models. 582 Journal of the Royal Statistical Society: Series B (Statistical Methodology) 583 63, 425464. doi: 10.1111/1467-9868.00294 . 584 Lin, Z.L., Beck, M.B., 2012. Accounting for structural error and uncertainty 585 in a model: An approach based on model parameters as stochastic pro586 cesses. Environmental Modelling & Software 27-28, 97111. doi: 10.1016/ 587 j.envsoft.2011.08.015 . 588 Liu, C., Zachara, J.M., 2001. Uncertainties of Monod kinetic parameters 589 nonlinearly estimated from batch experiments. Environmental Science & 590 Technology 35, 133141. doi: 10.1021/es001261b . 591 24 Ma²i¢, A., Srinivasan, S., Billeter, J., Bonvin, D., Villez, K., 2017. Shape 592 constrained splines as transparent black-box models for bioprocess mod593 eling. Computers & Chemical Engineering 99, 96105. doi: 10.1016/j. 594 compchemeng.2016.12.017 . 595 Murphy, K.P., 2012. Machine learning: A probabilistic perspective. MIT 596 Press, Cambridge, MA, USA. 597 Neumann, M.B., Gujer, W., 2008. Underestimation of uncertainty in sta598 tistical regression of environmental models: Inuence of model struc599 ture uncertainty. Environmental Science & Technology 42, 40374043. 600 doi: 10.1021/es702397q . 601 Olsson, G., 2012. ICA and me - A subjective review. Water Research 46, 602 15851624. 603 Parsons, S., 2001. Qualitative methods for reasoning under uncertainty. MIT 604 Press, Cambridge, MA, USA. 605 Petersen, B., Gernaey, K., Devisscher, M., Dochain, D., Vanrolleghem, P.A., 606 2003. A simplied method to assess structurally identiable parameters 607 in Monod-based activated sludge models. Water Research 37, 28932904. 608 doi: 10.1016/S0043-1354(03)00114-3 . 609 Rao, K.D., Kushwaha, H.S., Verma, A.K., Srividya, A., 2008. Epistemic 610 uncertainty propagation in reliability assessment of complex systems. In611 ternational Journal of Performability Engineering 4, 7184. 612 Reichert, P., Mieleitner, J., 2009. Analyzing input and structural uncertainty 613 of nonlinear dynamic models with stochastic, time-dependent parameters. 614 Water Resources Research 45, W10402. doi: 10.1029/2009WR007814 . 615 Reichert, P., Schuwirth, N., 2012. Linking statistical bias description to 616 multiobjective model calibration. Water Resources Research 48, W09543. 617 doi: 10.1029/2011WR011391 . 618 Renard, B., Kavetski, D., Kuczera, G., Thyer, M., Franks, S.W., 2010. Un619 derstanding predictive uncertainty in hydrologic modeling: The challenge 620 of identifying input and structural errors. Water Resources Research 46. 621 doi: 10.1029/2009WR008328 . 622 25