In engineering practice, mo del-based design requires not only a go o d pro cess- based mo del, but also a go o d description of sto c hastic disturbances and measuremen t errors to learn credible parameter v alues from observ ations. Ho w ev er, t ypical metho ds use Gaussian error mo dels, whic h often cannot de- scrib e the complex temp oral patterns of residuals. Consequen tly , this results in o v ercon dence in the iden ti ed parameters and, in turn, optimistic reactor designs. In this w ork, w e assess the strengths and w eaknesses of a metho d to statistically describ e these patterns with auto correlated error mo dels. This metho d pro duces increased widths of the credible prediction in terv als fol- lo wing the inclusion of the bias term, in turn leading to more conserv ativ e design c hoices. Ho w ev er, w e also sho w that the augmen ted error mo del is not a univ ersal to ol, as its application cannot guaran tee the desired reliabilit y of the resulting w astew ater reactor design. Keywor ds: bias description; kinetic mo del; pro cess design; w astew ater treatmen t; uncertain t y Email addr ess: villezk@ornl.gov (Kris Villez) Pr eprint submitte d to R eliability Engine ering and System Safety July 21, 2020 Cop yrigh t notice This man uscript has b een authored in part b y UT-Battelle, LLC, under con tract DE-A C05-00OR22725 with the US Departmen t of Energy (DOE). The US go v ernmen t retains and the publisher, b y accepting the article for publication, ac kno wledges that the US go v ernmen t retains a nonexclusiv e, paid-up, irrev o cable, w orldwide license to publish or repro duce the pub- lished form of this man uscript, or allo w others to do so, for US go v ern- men t purp oses. DOE will pro vide public access to these results of feder- ally sp onsored researc h in accordance with the DOE Public A ccess Plan (h ttp://energy .go v/do wnloads/do epublic- access-plan). 1. In tro duction 1 In curren t en vironmen tal engineering practice, deterministic pro cess-based 2 mo deling is a common to ol to b etter understand the functioning of complex 3 w astew ater collection and treatmen t systems. The gold standard is to im- 4 pro v e prediction p erformance of our mo dels b y tting them to observ ations. 5 Consequen tly , the adv en t of ubiquitous sensing leads to an unin tended y et 6 commonly observ ed situation where sensors rev eal more details than mec h- 7 anistic mo dels can capture. When this is the case, uncertain t y estimates 8 obtained from statistical inference with mec hanistic mo dels are almost cer- 9 tainly to o narro w as the applied mo del structure is to o restrictiv e relativ e 10 to the observ ed realit y . A long-standing question is therefore whether risk- 11 based design, based on uncertain t y estimates from statistical inference with 12 mec hanistic mo dels, is actually feasible. In this w ork, w e test one metho d 13 designed to address this issue to a case of WWTP design and discuss its 14 p oten tial and limitations. 15 A ccoun ting for mo del parameter uncertain t y is crucial for risk-based decision- 16 making, including infrastructure design and op erations (e.g., Cagno et al., 17 2011; Sc heidegger et al., 2013; Kabir et al., 2015; Sc heidegger et al., 2015; 18 Jensen and Jerez, 2018). Con v en tional metho ds for uncertain t y analysis are 19 based on a t w o-step approac h, consisting of (a) quan ti cation of input un- 20 certain t y , measuremen t uncertain t y , and subsequen t uncertain t y of mo del 21 parameters follo w ed b y (b) propagation of the quan ti ed uncertain t y to the 22 system p erformance measure of in terest (e.g., V an Griensv en and Meixner, 23 2007; Sin et al., 2009; Guo and Murph y, 2012; Del Giudice et al., 2016). 24 2 Ho w ev er, it has b een demonstrated b efore ho w systematic de ciencies in 25 mo del structure, next to input and measuremen t uncertain t y , also lead to bi- 26 ased mo del parameters and, consequen tly , incorrect design of infrastructural 27 elemen ts, suc h as biological reactor systems (Neumann and Gujer, 2008). 28 Unfortunately , to our kno wledge, no one has attempted to pro vide a metho d 29 to solv e this particular problem, i.e. to iden tify systematic discrepancies in 30 pro cess-based mo dels so to accoun t for them during mo del-based design. 31 Recen tly , statisticians ha v e b een suggesting a promising approac h to solv e 32 this dilemma. The underlying idea is to not assume iden tically and indep en- 33 den tly distributed (i.i.d.) errors for mismatc hes b et w een mo dels and obser- 34 v ations (Liu and Zac hara, 2001), but to explicitly accoun t for mismatc hes b y 35 adding a sto c hastic auto-correlated pro cess to the i.i.d. measuremen t error 36 mo del (Craig et al., 2001; Kennedy and O'Hagan, 2001; Ba y arri et al., 2007). 37 This is kno wn as the bias description metho d and fo cuses on the mo deling 38 of the symptoms of a mismatc h b et w een mo del structure and realit y . While 39 this do es not iden tify or tac kle the ro ot cause of these symptoms, it has 40 b een pro v en to b e a computationally e cien t to ol to increase the reliabilit y 41 of mo del-based predictions compared to standard regression approac hes in a 42 v ariet y of systems from lak es to natural catc hmen ts to urban h ydrology (Di- 43 etzel and Reic hert, 2012; Reic hert and Sc h u wirth, 2012; Del Giudice et al., 44 2015). Therefore, w e exp ect that the bias description metho d also impro v es 45 the reliabilit y of predictions with structurally de cien t w astew ater treatmen t 46 mo dels in view of risk-based design. Sp eci cally , adding a sto c hastic mea- 47 suremen t error term to a mo del giv en the same amoun t of exp erimen tal 48 measuremen ts is exp ected to reduce the relativ e information-ric hness of the 49 exp erimen tal data and lead to larger credibilit y in terv als of mo del parame- 50 ters, wider prediction in terv als and, b y a v oiding o v ercon den t predictions, a 51 more trust w orth y design. 52 Note that the bias description metho d can b e regarded as a grey b o x or 53 h ybrid mo delling strategy . Indeed, the resulting mo del consists of a mec h- 54 anistic mo del for the studied pro cess (white b o x) and a sto c hastic mo del 55 for auto-correlated measuremen t errors (blac k b o x). Other grey b o x ap- 56 proac hes ma y b e based on the inclusion of time-v arian t parameters (Reic hert 57 and Mieleitner, 2009; Lin and Bec k, 2012) or in tegration of non-parametric 58 elemen ts in to a mo del structure that is mec hanistic otherwise (Ma²i¢ et al., 59 2017). 60 In this con tribution, w e apply the bias description metho d to in v estigate 61 the impact of mo del structure de cits for pro cess design. W e use Neumann 62 3 and Gujer (2008) as a b enc hmark to ev aluate the b ene ts and limitations 63 of the bias description metho d for mo del-based design and refer to it as the 64 r efer enc e study . While this reference study concerns a conceptually simple 65 case, using it in this study highligh ts (a) that the apparen t simplicit y of this 66 case is rather deceptiv e and (b) that c hallenges asso ciated with mo del-realit y 67 mismatc h are to b e exp ected for b oth simple and complex systems. 68 2. Material and metho ds 69 2.1. Applie d err or mo dels 70 In a v ast ma jorit y of en vironmen tal mo deling studies, the measuremen t 71 error is assumed to b e i.i.d. F or example, Hauduc et al. (2015) compares an 72 extensiv e list of mo del p erformance criteria for w astew ater treatmen t mo d- 73 elling y et do es not list an y criterion whic h accoun ts for auto correlated mo del 74 prediction errors. One approac h considered in Cierk ens et al. (2012), consists 75 of do wnsampling time series to a v oid the app earance of auto correlation. As 76 explained in the same study , this leads to an ine cien t use of the a v ailable 77 data and, more imp ortan tly , cannot accoun t at all for mo del structure de cits 78 as a p oten tial ro ot cause of auto correlated residuals. Ignoring the presence 79 of auto correlated residuals w as sho wn to lead to o v ercon dence in the pro- 80 duced mo del and, subsequen tly , p o or decision-making, as w as sho wn also in 81 the reference study . Most often, a Gaussian distribution is assumed for the 82 measuremen t errors. Suc h a mo del of measuremen t error cannot accoun t for 83 systematic deviations b et w een the assumed mo del and the observ ed mea- 84 suremen ts, i.e. bias. One w a y of accoun ting for bias is b y adding terms, suc h 85 as a sto c hastic auto correlated error bias term, to the measuremen t equation 86 (Craig et al., 2001; Kennedy and O'Hagan, 2001; Ba y arri et al., 2007). In 87 this w ork, w e describ e the observ able output time-series (i.e., measured con- 88 cen tration, yo ) as a sum of a deterministic dynamic mo del output ( y , the 89 mo deled concen tration), a classical Gaussian measuremen t 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 mo del and 92 are those of the sto c hastic parts (errors). The bias term b( ) decrib es an 93 auto correlated error ( b( )N(0;b(b; )) ) and can b e included to accoun t 94 4 for time-dep enden t deviations b et w een mo del and observ ations (Reic hert and 95 Sc h u wirth, 2012). Note that this bias term represen ts a sto c hastic pro cess, 96 th us describing aleatory uncertain t y , although the deviations b et w een mo del 97 and observ ations ma y actually b e systematic, p ossibly ev en deterministic. 98 These deviations are exp ected to b e systematic when they are caused b y a 99 lac k of kno wledge ab out the true data-generating pro cess. This lac k of kno wl- 100 edge is t ypically c haracterized as a source of epistemic uncertain t y rather than 101 aleatory uncertain t y . 102 The bias term has t w o parameters, the standard deviation b and the 103 correlation length : 104 b(i;j) :=b2ejtitjj2= : (2) The random measuremen t error is temp orally indep enden t ( eN(0;e(e)) ) 105 and is c haracterized b y the parameter e : 106 e(i;j) :=( e2; i=j 0; i6=j: (3) T ogether, these error terms with parameters =f ;b;eg accoun t for 107 the fact that the deterministic mo del ma y not repro duce the mo deled data 108 set exactly . Note that the sym b ols b ande are c hosen to con v ey the idea 109 that they b oth describ e the magnitude of v ariation of a sto c hastic term in the 110 measuremen t equation. The sym b ol for the correlation length, , is c hosen 111 to highligh t the fact that it describ es a time-scale. 112 The statistical form ulation in Eq. 1 naturally leads to the lik eliho o d 113 functionL(yoj; ) whic h describ es ho w lik ely the considered mo del with 114 parameters (; ) generated the recorded data, yo . The lik eliho o d of the 115 measuremen ts conditional to the mo del parameters is: 116 L(yoj; ) =(2)n 2q det exp 1 2h yoyiT ()1h yoyi (4) where  is the v ariance-co v ariance matrix for the sto c hastic deviations 117 b et w een mo del and observ ations: 118 5  := e+ b; (5) with b ande de ned in 2 and 3, and with n equal to the n um b er of 119 measuremen ts. In order to in terpret the results of parameter estimation, w e 120 de ne the parameters and suc h that2 e:= (1 )2and2 b:= 2. 121 This means w e can express the v ariance-co v ariance matrix ab o v e equiv alen tly 122 as: 123  :=2h (1 )Inn+ Ki (6) K(i;j) :=ejtitjj2= (7) In this form,  is a measure for the o v erall spread of the deviations b e- 124 t w een the mo del and the measuremen ts and is a parameter that de nes 125 the relativ e imp ortance of the bias in the o v erall v ariance-co v ariance matrix. 126 Meaningful v alues for are b et w een 0 and 1, with = 1 leading to the omis- 127 sion of the indep enden t measuremen t noise ( e= 0 ) and = 0 expressing 128 that there is no bias ( b= 0 ). Note that setting = 0 repro duces the mo del 129 without a bias term. Put otherwise, the mo del with bias term includes the 130 mo del without bias term as sp ecial case. 131 The addition of a bias term accoun ts for underestimation of parameter 132 uncertain t y when a con v en tional y et unrealistic distribution for the mo del 133 error is assumed (e.g., uncorrelated). This is exp ected to pro duce a wider 134 predictiv e distribution, p ossibly leading to a b etter quan ti cation of and a 135 reduction of the risk of under-design or o v er-design. Ho w ev er, sp ecial atten- 136 tion m ust b e giv en to the parameter estimation metho d as increased mo del 137 exibilit y can lead to uniden ti abilit y (see e.g., Renard et al., 2010). 138 2.1.1. Mo del p ar ameter estimation 139 W e apply a Ba y esian approac h for t w o reasons. First, w e fa v our a Ba y esian 140 framew ork as a w a y to mak e prior b eliefs explicit. Second, without an y form 141 of prior, some of the parameters of the v ariance-co v ariance matrix  can b e 142 structurally uniden ti able (for de nitions, see Do c hain et al., 1995; Do c hain 143 and V anrolleghem, 2001; P etersen et al., 2003). More sp eci cally , when = 0 144 the matrix K equals the iden tit y matrix and lik eliho o d L(yoj; ) b ecomes 145 insensitiv e to the v alue of . As a result, no unique v alue for can b e 146 6 iden ti ed under an y circumstances as long as = 0 , i.e. is structurally 147 uniden ti able. In the form ulation with e andb , an y increase of e can 148 b e comp ensated exactly b y an equiv alen t decrease of b when = 0 . F or 149 small v alues of , e.g. close to the measuremen t in terv al or smaller, this is 150 exp ected to lead to a lac k of practical iden ti abilit y , ev en if structural iden ti- 151 abilit y could b e guaran teed in principle. In early exp erimen ts with uniform 152 priors for , w e observ ed that this can induce a lac k of con v ergence and p o or 153 mixing conditions for the applied sampling metho ds, similar to observ ations 154 describ ed in Renard et al. (2010). Applying an informativ e prior solv es this 155 iden ti abilit y problem and can therefore also b e in terpreted as a form of 156 regularization (e.g., Scales and T enorio, 2001; Murph y, 2012; Hastie et al., 157 2015). 158 Ba y esian calibration aims at c haracterizing the distribution describ ed b y 159 the p osterior lik eliho o d L(; jyo)/ L(yoj; )L(; ) , where the prior 160 lik eliho o dL(; ) expresses the prior b eliefs ab out the parameters. In this 161 w ork, the p osterior distribution is appro ximated with a sample of L(; jyo) 162 dra wn with a Mark o v Chain Mon te Carlo (MCMC) sampler (see Numerical 163 implemen tation). 164 2.2. Biokinetic mo del p ar ameter identi c ation with b atch exp eriments 165 T o study the eects of mo del structure error and the utilit y of the bias de- 166 scription metho d, w e execute sim ulations with the dynamic biokinetic mo del 167 used in the reference study . Concretely , a series of batc h exp erimen ts is sim- 168 ulated in whic h a substrate, with concen tration s(t) , is consumed b y a cell 169 culture with a xed concen tration. The con v ersion rate r(t) dep ends on the 170 substrate b y means of time-in v arian t T essier kinetics so that one can write: 171 ds(t) dt=r(t) (8) s(t= 0) =s0 (9) r(t) =rTessier max 1exp s(t) KTessier S x(t) (10) During the exp erimen t, noisy measuremen ts of the true substrate con- 172 cen tration, yo(t) , are sim ulated b y the follo wing measuremen t error mo del, 173 whic h 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 eac h sim ulation are the same as in the reference 175 study:s0 (initial substrate concen tration, 5g=m3),rmax (maxim um con- 176 v ersion rate, 1g=m3:h ). The sim ulated time is T= 8 hours. The a nit y 177 constan t (KTessier S ) and measuremen t error standard deviation ( e ) are v aried 178 y et constan t in ev ery sim ulated exp erimen t. KTessier S is v aried from 0:1g=m3179 to1:5g=m3in steps of 0:2g=m3. This allo ws sim ulating a wide range of 180 pro cess conditions, including b oth lo w and high v alues for KTessier S relativ e 181 to the initial substrate concen tration. T w o v alues for the sim ulated e are 182 considered, as in the reference study . In the lo w-noise case, e tak es the v alue 183 0:01g=m3. In the high-noise case, it tak es the v alue 0:1g=m3. The v ector 184  equals s0;rmax;KTessier ST. The t w o sim ulated noisy time series obtained 185 withKTessier S = 0:7g=m3are sho wn in the supplemen tary information (Fig. 186 S.1 and Fig. S.2). 187 F or eac h sim ulation exp erimen t, parameter iden ti cation is executed with 188 four distinct mo del structures. The rst mo del matc hes the ab o v e mo del 189 structure (Eq. 8-Eq. 12) exactly . This represen ts an idealized situation where 190 the structure of the calibrated mo del matc hes realit y (ground truth) exactly . 191 The iden ti ed parameters are S0 ,Tessier max ,KTessier S , ande . A second mo del 192 is obtained b y replacing the T essier kinetics with the alternativ e and more 193 commonly used Mono d kinetics. Practically , Eq. 10, is replaced with the 194 follo wing equation: 195 r(t) =rMonod maxs(t) KMonod s +s(t)(13) The estimated parameters are no w s0 ,Monod max ,KMonod s , ande with 196 = s0;rmax;KMonod ST. This case represen ts the lik ely situation that a 197 mo deling practitioner uses the common-place Mono d mo del structure and 198 do es not observ e the bias that results. This is v ery lik ely in the high-noise 199 case (see reference study). Giv en this di cult y , the sto c hastic bias term de- 200 scrib ed ab o v e is included to capture the systematic deviations b et w een the 201 mo del predictions and measuremen ts. T o ac hiev e this, the previously applied 202 measuremen t equation (Eq. 11) is replaced with the follo wing 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 ab o v e and b ande reparametrized with 204 and . This results in a third mo del, where a T essier mo del is com bined 205 with the statistical bias description and whic h requires sp eci cation of the 206 parameters s0 ,Tessier max ,KTessier S , , , and . As the T essier mo del has 207 the same structure as the data-generating mo del, one can exp ect a go o d 208 mo del t with close to 0 and estimates of s0 ,Tessier max ,KTessier S , and 209 that are close to ground truth v alues. Finally , the fourth mo del com bines 210 the presumed Mono d kinetics with the statistical bias description and the 211 iden ti ed parameters are s0 ,Monod max ,KMonod s , , , and . In this case w e 212 can exp ect that the presen t mo del structure bias is accommo dated b y means 213 of the statistical bias description. If so, this should increase the width of the 214 prediction in terv als and thereb y impro v e the reliabilit y of the mo del (Reic hert 215 and Sc h u wirth, 2012). 216 2.3. Numeric al implementation 217 The biokinetic mo del, parameter estimation, and uncertain t y propaga- 218 tion w ere implemen ted in Matlab (R2019a). The prior probabilities for the 219 parameters w ere set based on the authors' exp erience. They are all indep en- 220 den t of eac h other. All priors are uniform, except for  and . The prior 221 lik eliho o d for  is prop ortional to its in v erse and is equiv alen t to the Jereys 222 prior conditional to xed v alues for all other parameters (see Bo x and Tiao, 223 1973). The prior lik eliho o d for is the sine function supp orted b et w een 0 and 224 2T . This prior equals zero at = 0 and = 2T and one at =T . This ex- 225 presses the sub jectiv e b elief that the auto correlation length of the deviations 226 due to mo del structure error is exp ected to b e similar to the duration of the 227 exp erimen t. The priors are sp eci ed completely in T able 1. W e rst run an 228 adaptiv e MCMC algorithm (Vihola, 2012) to nd a go o d guess for the maxi- 229 m um a p osteriori estimates and a go o d prop osal v ariance-co v ariance matrix. 230 With these results, w e execute a (non-adaptiv e) MCMC algorithm to obtain 231 20;000 samples fromL(; jyo) . The rst 10;000 samples are considered to 232 corresp ond to the burn-in phase of the sampler, during whic h eects of the 233 initial sample ma y still b e apparen t. These samples are therefore discarded, 234 as is common in practice (Gilks et al., 1996). 235 9 T able 1: Prior probabilit y distributions for the parameters. P arameter Distribution Lo w er b ound Upp er b ound Unit s0 Uniform 0 +1g=m3 rmax Uniform 0 +1g=m3:h KS Uniform 0 +1g=m3  In v erse (/1 ) e12e+12g=m3 Uniform 0 1  Sine (/sin( 2 T) ) 0 2 T h 2.4. Using information gaine d in p ar ameter estimation for design 236 The obtained parameter estimates for eac h mo del are insp ected b y means 237 of visual insp ection. In addition, w e test the reliabilit y of the calibrated 238 mo dels for subsequen t pro cess design. T o this end, the steady state substrate 239 concen tration is computed for a CSTR with cells gro wing according to the 240 aforemen tioned ground truth T essier kinetics. The dilution rate D is set at 241 0:5h1, corresp onding to a h ydraulic residence time of 2 h. The ground 242 truth steady state substrate concen tration can b e computed as follo ws: 243 s(t!1 )CSTR =KTessier slog 1D rmax (17) The ab o v e equation is also used to sim ulate the mo del-based concen tra- 244 tion estimate b y replacing the true v alue of KTessier s with its estimates. T o 245 do the same with the mo dels exhibiting Mono d kinetics, the steady state 246 concen tration estimate is computed as follo ws: 247 s(t!1 )CSTR =DKMonod s rmaxD(18) The steady state concen trations are computed assuming that a p erfectly 248 accurate v alue for rmax is a v ailable from a separate w ashout exp erimen t (as in 249 the reference study). The a nit y constan t, KMonod s orKTessier s , is therefore 250 the only parameter whose v alues are based on the batc h exp erimen t. In all 251 cases, the reliabilit y of the predicted steady-state concen tration in the CSTR 252 is assessed b y visual insp ection. Considering that the p erfect information 253 ab outrmax is not accoun ted for during estimation of the a nit y constan t, w e 254 10 also pro vide results in the supplemen tary information for a mo di ed v ersion 255 of our metho d. In this case, w e assume no w ashout exp erimen t has b een 256 executed and use the sampled v alues for rmax in Eq. 17 and Eq. 18 instead 257 of the assumed ground truth v alue. 258 2.5. Softwar e availability 259 All data and n umerical metho ds to repro duce our results are a v ailable as a 260 stand-alone pac k age for the Matlab (R2019b) platform. The v ersion used for 261 this article is added to the supplemen tary materials. The most recen t v ersion 262 can b e found on Gitlab ( h ttps://gitlab.com/krisvillez/biasdescription ). 263 3. Results 264 3.1. Par ameter estimation 265 The eect of bias and the use of the statistical bias description are demon- 266 strated rst b y means of the marginal p osterior parameter distributions ob- 267 tained with data sim ulated with the T essier mo del with KTessier s = 0:7g=m3268 ande= 0:01g=m3. In Fig. 1 the empirical cum ulativ e densit y functions 269 are sho wn for the idealized case where the mo del structure is correct. The 270 results without bias description are sho wn with a dashed red line. One can 271 observ e easily that the estimates of the deterministic part of the mo del are 272 fairly accurate. Indeed, the distributions are b oth narro w (all relativ e stan- 273 dard deviations are b elo w 2% in magnitude) and close to the true v alues 274 (all median estimates are within 1:5% of the ground truth). This is not 275 surprising. In con trast, the median estimate of  is ab out 35% smaller than 276 the sim ulated v alue while its precision remains small also (relativ e standard 277 deviation: 1.1%). 278 The red dashed lines in Fig. 2 sho w the empirical cum ulativ e densit y 279 functions pro duced b y using the Mono d mo del structure while using the 280 con v en tional error mo del without bias description. This results in biased 281 estimates of the con v ersion rate ( rMonod max ). These estimates do not re ect 282 the v alues used in the ground truth sim ulation, as discussed in the reference 283 study . F or instance, the relativ e deviation median estimates are +19% for 284 rMonod max . This is explained as a consequence of using a dieren t rate function 285 during estimation. Indeed, the parameter v alues for rmax (andKS ) are used 286 to comp ensate for the mo del structure error. The relativ e standard devia- 287 tions for all parameters remain b elo w 1.1% in magnitude, except for KMonod s288 (2.7%). 289 11 Figure 1: Batc h exp erimen t sim ulated with the T essier mo del and Ks= 0:7g=m3and e= 0:01g=m3. Marginal p osterior distributions for the calibrated T essier mo del without bias description (dashed red line) and with bias description (full blue line). V ertical full lines indicate the true data-generating parameter v alues ( rmax ,KS ,s0 , , ). Blac k dotted lines indicate the informativ e prior ( , ). This is dieren t in the case where the Mono d mo del structure is used 290 for estimation (Fig. 2). As exp ected, the p osterior distribution for is no w 291 lo cated to the righ t of its prior and indicates a large con tribution of bias 292 to the prediction error. A t the same time, the p osterior for is to the 293 left of the prior, whic h means the auto correlation length is shorter than 294 the exp erimen tal time length (8h). The residuals obtained with maxim um 295 lik eliho o d estimation of the Mono d mo del no w app ear correlated as w ell (see 296 Fig. S.2 in the supplemen tary information). 297 12 Figure 2: Batc h exp erimen t sim ulated with the Mono d mo del and Ks= 0:7g=m3and e= 0:01g=m3. Marginal p osterior distributions for the calibrated Mono d mo del without bias description (dashed red line) and with bias description (full blue line). V ertical full lines indicate the data-generating parameter v alues in the T essier mo del used for sim ulation (rmax ,KS ,s0 , , ). Blac k dotted lines indicate the informativ e prior ( , ). 3.2. Mo del-b ase d pr e diction: use of p ar ameter estimates for CSTR design 298 F or the purp ose of mo del-based design it is imp ortan t to predict the 299 steady state concen tration prop erly for a giv en desired dilution rate. As 300 indicated ab o v e, the steady state concen tration is predicted at a dilution 301 rate of 0:5h1for b oth the ground truth as w ell as with the four constructed 302 mo dels. 303 The lo w-noise case is discussed rst. The iden ti ed Mono d mo del pa- 304 rameter sets are used to predict the steady-state concen tration b y means of 305 Eq. 18. This is rep eated for ev ery sim ulated v alue of KTessier S . The ratios of 306 the predicted steady-state concen trations to the true steady-state concen tra- 307 13 tion are visualized in Fig. 3. It is clear that using the Mono d mo del leads 308 to a signi can t bias in these predictions. This situation is ho w ev er easy to 309 iden tify b y insp ection of the p osterior of (as explained ab o v e). In addi- 310 tion, mo del structure error in the lo w-noise case is also iden ti ed easily b y 311 means of frequen tist metho ds (cfr. reference study). W e therefore assume 312 that the mo deler applies one of these to ols and thereb y successfully iden ti es 313 the presence of bias. 314 Figure 3: Distributions of the ratio of predicted steady state concen trations to the ground truth concen tration as a function of the a nit y constan t ( KS )  Lo w noise case ( e= 0:01g=m3). Red horizon tal whisk ers indicate the t w o-sided 99% credible in terv als. Results are sho wn for the Mono d mo del without bias description. As the ideal ratio of 1 falls outside the credibilit y in terv als, w e see that using a wrong mo del for reactor design will lead to a failure to ac hiev e the desired p ollutan t remo v al capacit y . This situation is easily detected in this lo w-noise case, also when using traditional metho ds. As the high-noise case is more c hallenging for frequen tist metho ds, w e 315 14 discuss it in greater detail. Fig. 4 sho ws all results obtained for ev ery assumed 316 structure for the deterministic mo del (Mono d/T essier) and for the sto c hastic 317 mo del parts (with and without bias). The top panel displa ys the results 318 obtained with the T essier mo del. It is easy to see that the predictions are 319 un biased as w ell as precise when no bias description is incorp orated in the 320 error mo del (left-side b ean plots), except when KS= 0:1g=m3. W e return 321 to the latter case b elo w. If the sto c hastic bias term is added (righ t-side 322 b ean plots), the distribution of the predictions b ecomes wider without an y 323 meaningful shift of the a v erage predictions ( KS= 0:3g=m3and higher). 324 Consider that the pro cess design is based on the sho wn 99% credible limits. 325 As adding bias increases the predicted upp er limit, this will result in a larger 326 design v olume to accoun t for the p erceiv ed increase in uncertain t y . As suc h, 327 accoun ting for uncertain t y in the design w ould lead to a more conserv ativ e 328 design compared to the case without bias description. 329 The predictions obtained with a Mono d mo del are sho wn in the b ot- 330 tom panel of Fig. 4. Dep ending on the ground truth v alue for KTessier s , 331 rather sev ere under-prediction or o v er-prediction results when no bias term 332 is added (left-side b ean plots). The predictiv e uncertain t y (spread) is similar 333 to the case with the T essier mo dels without bias description. As a result, 334 the t w o-sided 99% credible in terv als include the ground truth in only one 335 of the sim ulated cases ( Ks= 0:5g=m3), indicating the lac k of reliabilit y 336 of the mo del-based prediction in terv als. The addition of a bias description 337 term abates this issue to some exten t. Ev en though bias is still presen t, the 338 addition of a phenomenological bias term leads to increased widths of the 339 credible in terv als and a more conserv ativ e design. Ho w ev er, the ground truth 340 is included in the 99% credible in terv als in t w o cases only ( Ks= 0:5g=m3, 341 Ks= 0:7g=m3). Th us, the bias description cannot ensure a reliable design 342 without further mo di cation. Insp ecting the p osterior of is more useful. 343 As one can see in Fig. 5, in all cases without mo del structure error (T essier 344 mo del), the p osterior of is shifted to w ards the left of its prior, th us sug- 345 gesting that the mo del structure is de ned w ell. F or the case with mo del 346 structure error (Mono d mo del), the p osterior is shifted to the righ t or re- 347 mains relativ ely close to the prior, except for KS= 0:1g=m3. One reason 348 is that this v alue is equal to the measuremen t error standard deviation. A 349 second lik ely factor is that the length of time during whic h the pro duced 350 exp erimen tal data are sensitiv e to the v alue of KS is v ery short. Indeed, 351 the substrate concen tration in the sim ulated exp erimen t is b et w een half and 352 t wice the v alue for KS for less than 3% of the duration of the exp erimen t. 353 15 Figure 4: Distributions of the ratio of predicted steady state concen trations to the ground truth concen tration as a function of the a nit y constan t ( KS )  High noise case ( e= 0:1g=m3). Red horizon tal whisk ers indicate the t w o-sided 99% credible in terv als. Left side b eans: without bias description; Righ t side b eans: with bias description. T op: T essier mo del - All credible in terv als include the ideal ratio (equal to 1), except for the sim ulation with KS= 0:1g=m3without bias term. The uncertain t y increases when a bias term is added to the mo del. Bottom: Mono d mo del  Including the bias term in the mo del increases the reliabilit y of the credible in terv als. These in terv als include the ideal ratio for t w o cases ( KS= 0:5 and 0:7g=m3). F orKS= 0:3g=m3this already amoun ts to 8.6%. Since the mo del structure 354 error primarily relates to the curv ature of the con v ersion rate in this region, 355 it follo ws that mo del structure error will alw a ys b e di cult to detect when 356 this time fraction is lo w. 357 In the supplemen tary information, w e pro vide results obtained with the 358 mo di ed metho d. W e omit the information obtained during the w ashout ex- 359 p erimen t during prediction. In this case, the uncertain t y in the predictions 360 16 Figure 5: Distributions of the parameter for b oth mo dels with a bias term in all high- noise cases ( e= 0:1g=m3). When the T essier mo del is selected (no mo del structure error), the p osterior distribution of is shifted to the left of the prior, th us suggesting the kinetic mo del structure is adequate. In con trast, the p osterior of is similar to or lo cated at the righ t of the prior when the Mono d mo del is used in all but one case ( KS= 0:1g=m3), th us pro viding a useful indication of mo del structure error. is reduced signi can tly to the p oin t that none of the 99% credible in terv als 361 include the ground truth (see Fig.S.3). This is explained b y the fact that 362 the estimates for rmax andKS exhibit strong correlation (see supplemen- 363 tary information for details). Ho w ev er, since the mo di cation relates to the 364 prediction step only , one can still use the p osteriors for as a detection 365 mec hanism for bias. 366 17 4. Discussion 367 4.1. Summary and limitations of the exp erimental simulation study 368 Summary. The n umerical results describ ed ab o v e suggest that the inclusion 369 of an additiv e auto-correlated error pro cess in to a measuremen t error mo del 370 can impro v e the reliabilit y of mo del-based designs. This is true ev en when 371 only a subset of the iden ti ed parameters are used during prediction (here w e 372 only used the estimates for KS ) and ev en when the exp erimen tal setting for 373 prediction (steady state) is dieren t from the exp erimen tal conditions used 374 for mo del iden ti cation (batc h exp erimen t). In our case, the bias description 375 metho d impro v es the reliabilit y in all cases. Despite this impro v emen t, the 376 computed credible in terv als include the ground truth v alue only in a lim- 377 ited n um b er of cases with mo del structure error, meaning that guaran teed 378 reliabilit y cannot b e obtained with the studied metho d. Th us, the inclusion 379 of a bias description term for the purp ose of prediction can b e advised as 380 a relativ ely fast and easy w a y to accoun t for errors in the prop osed mo del 381 structure, ho w ev er only when one is unable to mo dify the mo del structure 382 itself. This is esp ecially relev an t in engineering applications where one is re- 383 stricted to sp eci c pro cess represen tations (e.g., Mono d kinetics) or soft w are 384 with limited exibilit y . While the bias description metho d impro v es the re- 385 liabilit y of the mo del predicitions only in a limited w a y , it is v ery useful as 386 a to ol to detect the presence of bias during mo del iden ti cation, esp ecially 387 when reform ulated with the parameter. 388 Limitations. In this study , a simple case w as c hosen delib erately for t w o 389 reasons. First, this enabled an ob jectiv e comparison of the bias descrip- 390 tiv e metho d with the historical results in the reference study (Neumann and 391 Gujer, 2008). Second, the apparen t simplicit y of the case also highligh ts the 392 c hallenge of generating reliable predictions with mec hanistic mo dels, induced 393 b y the t ypical lac k of exibilit y of suc h mo dels. The c hosen scop e also means 394 that our study comes with some limitations, whic h are: 395  The general applicabilit y of the bias description metho d is not demon- 396 strated. Ho w ev er, the bias description metho d could easily b e adapted 397 to more complex systems. One could incorp orate a bias term to ex- 398 press correlation b et w een m ultiple measuremen ts, of the same or dis- 399 tinct v ariables measured in the same lo cation or dieren t lo cations. In 400 this case, the co v ariance b et w een t w o measuremen ts, as expressed b y  , 401 w ould not only b e a function of (a) the time dierence ( titj , see (6)), 402 18 as in our study , but also of (b) spatial distance in one or more dimen- 403 sions and (c) eects of measuremen t error correlation b et w een distinct 404 sensors measuring the same or distinct v ariables. This generalization 405 of the presen t mo del is lik ely most con v enien t when the bias error term 406 is mo delled as spatio-temp oral Gaussian pro cess (e.g., De Cesare et al., 407 2001; Gneiting, 2002; Stein, 2005). 408  The metho ds applied in b oth the reference study and ours are based on 409 metho ds that accoun t for aleatory uncertain t y only . Ho w ev er, the lac k 410 of kno wledge ab out the mo del structure is t ypically epistemic in nature 411 and ma y therefore b e di cult to accoun t for in this w a y . Epistemic 412 uncertain t y ma y ho w ev er b e reduced b y using more exible mo dels 413 (Ma²i¢ et al., 2017) while increasing parametric uncertain t y , whic h can 414 b e handled as an aleatory source of uncertain t y with curren tly a v ailable 415 metho ds. Still, the adoption of alternativ e framew orks for uncertain t y 416 analysis (P arsons, 2001; Rao et al., 2008) ma y b e suited to handle 417 epistemic uncertain t y directly . In summary , the handling of epistemic 418 uncertain t y deserv es more atten tion. 419 4.2. Gener al c onse quenc es for pr actic al unc ertainty and r eliability analysis 420 Utility of the bias description metho d. In our opinion, the detection of sys- 421 tematic deviations b et w een the assumed mo del structure and the data-generating 422 pro cess is the most useful feature of the bias description metho d. F or this 423 reason, w e recommend that a mo del is insp ected for bias b y (a) adding a bias 424 term in the assumed mo del, sp ecifying a prior for alpha concen trated around 425 a strictly p ositiv e v alue, as suggested here, and (b) insp ecting the p osterior 426 of whenev er an inappropriate mo del structure is susp ected. Reform ulation 427 of the error mo del (bias + measuremen t error) with , , and pro v ed v ery 428 helpful as it enables in terpreting as an indicator for the relativ e imp or- 429 tance of mo del structure error. In cases where the p osterior probabilit y mass 430 is not shifted to w ards zero, relativ e to the prior, the mo deler should susp ect 431 the presence of bias. When this is detected, p oten tial mo del impro v emen ts 432 ma y include the use of time-dep enden t parameters (Reic hert and Mieleitner, 433 2009; Lin and Bec k, 2012) and/or input errors (Del Giudice et al., 2016) or 434 a c hange in mo del structure (Del Giudice et al., 2015; Ma²i¢ et al., 2017). 435 While the metho d increases the reliabilit y of the obtained steady-state p ol- 436 lutan t concen tration predictions, it is imp ortan t to note that the observ ation 437 of this b ene t dep ends strongly on the ro ot cause of the observ ed bias. F or 438 19 this reason, detection of bias should b e follo w ed b y exploratory analysis of 439 the residuals and dev elopmen t of a b etter mo del structure (e.g., Reic hert 440 and Mieleitner, 2009; Del Giudice et al., 2013). W e do not recommend ex- 441 ploiting the bias term for prediction without searc h for the underlying causes 442 for mo del structure de cits, esp ecially considering that the ground truth is 443 rarely included in the pro duced credible in terv als. Ultimately , the utilit y 444 of an y approac h dep ends on whether it can successfully describ e the rele- 445 v an t sources of the deviations b et w een mo del predictions and the measured 446 v ariables (Brynjarsdóttir and O'Hagan, 2014; W ani et al., 2019). 447 Par ameter interpr etation and tr ansfer ability. The mec hanistic in terpretation 448 of iden ti ed v alues for the parameters in the deterministic part of the mo del 449 is nearly imp ossible when bias is presen t. A dding an auto-correlated additiv e 450 error term con tributes to a b etter reliabilit y of the mo del predictions but can- 451 not pro vide a clearer in terpretation of the parameter v alues or a direction to 452 a more appropriate mo del structure. Indeed, the parameter estimates remain 453 biased. Imp ortan tly , this is a lik ely scenario in w astew ater engineering due to 454 the extremely simpli ed represen tation of biological pro cesses during mo del 455 construction. F urthermore, obtaining pro ofs of a lac k of bias is extremely 456 di cult to ac hiev e so that a straigh tforw ard in terpretation of parameter v al- 457 ues is unlik ely , ev en when the mo del structure ma y b e appropriate. Ho w ev er, 458 grey-b o x or h ybrid mo dels ma y oer in tepretabilit y and transparency at the 459 cost of computational eorts (see in tro duction ab o v e). 460 Data quality. The qualit y of the sim ulated measuremen ts in the studied case 461 is fairly high relativ e to curren t exp erience in the w astew ater sector. Ho w ev er, 462 sensor hardw are has b ecome increasingly robust in the last three decades 463 (Olsson, 2012) and there is no ob vious reason wh y this trend should stop 464 no w. It is therefore reasonable to exp ect that the presence of bias can b e 465 detected easily in the future, either b y statistical tests for auto-correlation 466 of the residuals, as in the reference study , or with descriptiv e metho ds, as in 467 this study . This will also facilitate the mo di cation of the mo del structure 468 in accordance to the en visioned high-qualit y data. 469 4.3. F utur e work 470 Through this w ork, w e iden ti ed sev eral a v en ues of further researc h. 471 These include: 472 20  Dev elop and study metho ds for parameter estimation and parameter 473 in terpretation under presence of mo del structure error. 474  Dev elop a systematic approac h to the form ulation of prior distributions, 475 esp ecially when exibilit y is at o dds with mo del structure or parameter 476 iden ti abilit y . 477  Ev aluation of exp erimen tal design metho ds to impro v e the c hances of 478 detection of mo del structure errors. 479  A dopt and ev aluate metho ds to handle epistemic uncertain t y in mo del- 480 based pro cess design and op eration. 481 5. Conclusions 482 In this pap er, w e in v estigated the c hallenge of structural mo del de cits in 483 risk-based reactor design. This is a relev an t problem, b ecause digitalization 484 will impro v e sensor resolution and spatial co v erage of reactors, whic h will 485 rev eal mismatc hes (i.e, bias) in our common engineering mo dels (whic h ha v e 486 b een dev elop ed in the data-scarce past, often b y grab sampling). Auto- 487 correlated mathematical form ulations ha v e b een suggested to impro v e the 488 description of suc h biases. 489 In summary , our study sho ws that 490  A dding auto-correlation terms in the measuremen t error mo del as a 491 w a y to accoun t for mo del structure de cits signi can tly impro v es the 492 reliabilit y of biokinetic mo dels. 493  Bias description enables accoun ting for predictiv e uncertain t y during 494 pro cess design to a large degree. This do es not pro duce a guaran teed 495 reliabilit y of the resulting design ho w ev er. It is therefore not a bullet- 496 pro of solution to the presence of mo del-realit y mismatc h. 497  The studied bias description metho d is an adequate to ol to iden tify the 498 presence of mo del structure de cits in presence of noisy exp erimen tal 499 data. 500 21 A c kno wledgmen ts 501 W e thank P eter Reic hert and Sanda Dejanic for their helpful insigh t in to 502 the studied problem. Marc B. Neumann ac kno wledges nancial supp ort pro- 503 vided b y the Spanish Go v ernmen t through the BC3 María de Maeztu ex- 504 cellence accreditation 2018-2022 (MDM-2017-0714) and the Ramón y Ca jal 505 gran t (R YC-2013-13628); and b y the Basque Go v ernmen t through the BER C 506 2018-2021 program. 507 6. Author con tributions 508 Dario del Giudice: Metho dology , Soft w are, Data analysis, W riting - Re- 509 view & Editing. Marc B. Neumann: Metho dology , Soft w are, Data analysis, 510 W riting - Review & Editing. Jörg Riec k ermann: Conceptualization, Metho d- 511 ology , W riting - Review & Editing. Kris Villez: Metho dology , Soft w are, 512 F ormal analysis, Data analysis, Visualization, W riting - Original Draft. 513