An integrated mechanistic modeling of a facultative pond: Parameter estimation and uncertainty analysis

Accepted Manuscript An integrated mechanistic modeling of a facultative pond: parameter estimation and uncertainty analysis Long T. Ho, Andres Alvarado, Josue Larriva, Cassia Pompeu, Peter Goethals PII:

S0043-1354(18)31035-2

DOI:

https://doi.org/10.1016/j.watres.2018.12.018

Reference:

WR 14313

To appear in:

Water Research

Received Date: 26 June 2018 Revised Date:

14 December 2018

Accepted Date: 14 December 2018

Please cite this article as: Ho, L.T., Alvarado, A., Larriva, J., Pompeu, C., Goethals, P., An integrated mechanistic modeling of a facultative pond: parameter estimation and uncertainty analysis, Water Research, https://doi.org/10.1016/j.watres.2018.12.018. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT

An integrated mechanistic modeling of a facultative pond: parameter

2

estimation and uncertainty analysis

3

Long T. Ho1*, Andres Alvarado2,3, Josue Larriva4,5, Cassia Pompeu1, Peter

4

Goethals1

RI PT

1

1

Department of Animal Sciences and Aquatic Ecology, Ghent University, Ghent, Belgium.

6

2

Departamento de Recursos Hídricos y Ciencias Ambientales, Universidad de Cuenca, Av. 12 de

7

Abril s/n, Cuenca, Ecuador.

8

3

9

4

SC

5

M AN U

Facultad de Ingeniería, Universidad de Cuenca, Av. 12 de Abril s/n, Cuenca, Ecuador. ETAPA, Empresa Pública Municipal de Telecomunicaciones, Agua Potable, Alcantarillado y

10

Saneamiento de Cuenca, Panamericana Norte km. 5 1/2, Ucubamba, Cuenca, Ecuador.

11

5

12

* Corresponding author. E-mail: [email protected]

AC C

EP

TE D

Facultad de Ciencia y Tecnología, Universidad del Azuay, Av. 24 de Mayo 7-77, Cuenca, Ecuador.

1

ACCEPTED MANUSCRIPT Abstract

14

Imitating natural lakes, pond treatment systems inherit a high complexity with interconnected web of

15

biochemical reactions and complex hydraulic processes. As such, its simulation requires a large and

16

integrated model, which has been a challenge for pond engineers. In this study, we develop an all-

17

encompassing model to gain a quantitative and comprehensive understanding of the hydraulic,

18

physicochemical and microbiological conversion processes in the most common pond, a facultative

19

pond. Moreover, to deal with an evitable issue of large mechanistic models as overparameterization

20

leading to poor identifiability, a systematic parameter estimation was implemented. The application of

21

sensitivity analysis reveals the most influential parameters on pond performance. Particularly, physical

22

parameters, such as vertical eddy diffusivity, water temperature, and maximum growth rate of

23

heterotrophs induce the most changes of organic matters while microbial assimilation and ammonia

24

volatilization appear to be main processes for nutrient removal. In contrast, the efficiency of phosphate

25

precipitation and nutrient biological removal via polyphosphate accumulating organisms and

26

denitrifying bacteria is limited. Identifiability problems are addressed mainly by the characterization

27

of light dependence of algal growth, interaction between water temperature and its coefficient, and the

28

growth of autotrophic bacteria while based on the determinant measures, the most important parameter

29

subsets affecting model outputs are related to physical processes and algal activity. After the

30

establishment of the influential and identifiable parameter subset, an automatic calibration with the

31

data collected from Ucubamba pond system (Ecuador) demonstrates the effect of high-altitude

32

climatic conditions on pond behaviors. An aerobic prevailing condition is observed as a result of high

33

light intensity causing accelerated algal activities, hence, leading to the limitation of hydrolysis,

34

anaerobic processes, and the growth of anoxic heterotrophs for denitrification. Furthermore, the output

35

of uncertainty analysis indicates that a large avoidable uncertainty as a result of vast complexity of the

36

applied model can be reduced greatly via a systematic approach for parameter estimation.

AC C

EP

TE D

M AN U

SC

RI PT

13

37

Keywords: Waste stabilization pond; Integrated mechanistic model; Sensitivity analysis;

38

Identifiability analysis; Parameter estimation; Uncertainty analysis;

2

ACCEPTED MANUSCRIPT 1. Introduction

40

Waste stabilization ponds (WSPs) are large shallow lagoons surrounded by earth embankments, where

41

wastewater is purified by completely natural processes. Taking advantage of algal photosynthetic

42

oxygenation, WSPs are a major choice for wastewater treatment plants (WWTPs) in many countries,

43

e.g. the US (>8,000 installations), France (>2,500 installations), and Germany (3,000 installations)

44

(EPA 2011). Despite being categorized as simple and basic operational system, processes inside the

45

pond systems are complicated. Facultative ponds (FPs), which are normally used as primary or

46

secondary treatment in pond series, display a high complexity because of a synergy between

47

microalgae and bacteria. This mutualistic relationship is one of the main reasons for the variation of

48

vertical oxygen profile, creating simultaneous occurrence of aerobic, anoxic, and anaerobic zones

49

inside FPs. To simulate this interconnected web of biochemical processes and reactions, an integrated

50

modeling approach is required, however, has been a challenge for pond engineers (Sah et al. 2011).

51

There are two main goals of applying models in environmental sciences in general. Firstly, models are

52

a powerful tool to describe and analyze the behavior of the systems. Second, models can be used to

53

forecast the responses of the systems under varying conditions. Each of these two goals addresses

54

different pathways of modeling techniques. For descriptive models, the task is to construct a simple

55

model with good estimates of model parameters which can simulate measured output with reasonable

56

accuracy (Omlin and Reichert 1999). A more complex model, which is unable to increase

57

considerably its goodness of fit or gain a better understanding of the system, is normally considered a

58

wasteful expense according to the parsimony principle of Spriet (1985). On the other hand, if the goal

59

of the model is to have a reliable prediction and reasonable approximation of model uncertainty,

60

simple model structure may result in poor predictive performance and significant underestimation of

61

uncertainty estimates (Reichert and Omlin 1997). From this perspective, an integrated mechanistic

62

model is a more suitable candidate. In fact, during the last two decades, there has been a shift from

63

data-driven models to mechanistic models due to the great interest in seeking better cost-effective

64

designs, process optimization, and pond upgrades (Ho et al. 2017). Besides, pond engineers also desire

65

to extrapolate their models to other systems under different driving conditions which is impossible

AC C

EP

TE D

M AN U

SC

RI PT

39

3

ACCEPTED MANUSCRIPT without simulating the causal mechanisms of the systems (Brun et al. 2001). However, pond engineers

67

are currently facing two fundamental problems in simulating this complex system. First, there is no

68

standard model which includes hydrodynamics, physicochemical as well as microbial conversion

69

processes (Sah et al. 2012). More importantly, these large and complex models are almost always

70

overparameterized, leading to a problem of poor identifiability as the number of available data for

71

parameter estimation is of significantly lower order than model complexity (Omlin and Reichert

72

1999).

73

Identifiability problem relating to the application of large mechanistic models has long been studied

74

and widely reported by the scientific community of WWTPs (Brun et al. 2002, Van Veldhuizen et al.

75

1999, Weijers and Vanrolleghem 1997). However, it is striking that the number of research paper

76

dealing with this issue in WSPs remains extremely low (Ho et al. 2017). In fact, to the author’s

77

knowledge, there have been only three mechanistic models in pond modeling, i.e. Kayombo et al.

78

(2000), Dochain et al. (2003), Beran and Kargi (2005) that conducted model calibration, but without

79

the consideration of the identifiability problem. In this study, we develop an all-encompassing

80

mathematical model which simulates not only biogeochemical processes of carbon, nitrogen, and

81

phosphorus removal but also the hydraulic and physical processes of a FP, the most common ponds in

82

WSP systems. We explicitly consider the interplay between biogeochemical and hydraulic submodels

83

in the mass balance of multiple variables. The main purpose is to gain a quantitative and

84

comprehensive understanding of the hydraulic, physicochemical and microbiological conversion

85

processes in the FP. This leads to the fact that many relevant processes, which are not yet included in

86

the previous models (Colomer and Rico 1993, Gehring et al. 2010, Sah et al. 2011), such as

87

phosphorus removal processes by phosphorus-accumulating organisms (PAO) or precipitation,

88

ammonia volatilization, sedimentation and resuspension, are included in the model. Subsequently, to

89

deal with the inevitable problem of overparameterization in this large mechanistic model, a systematic

90

procedure of parameter estimation was conducted. Firstly, sensitivity (SA) and identifiability analyses

91

(IA) are implemented via a practical technique proposed by Brun et al. (2002) where two criteria of

92

model parameters are addressed. First, the model output needs to be sufficiently sensitive to the

AC C

EP

TE D

M AN U

SC

RI PT

66

4

ACCEPTED MANUSCRIPT changes of the parameters. Second, these changes may not be compensated by the changes of other

94

parameters. As a result, a subset of influential and identifiable parameters is revealed, which,

95

subsequently, is calibrated by fitting the model to the data collected in Ucubamba WSP, Cuenca,

96

Ecuador in 2013. Finally, Monte Carlo simulation, a nonlinear error propagation with Latin

97

Hypercube Sampling technique, is applied for uncertainty analysis.

98

2. Methods

99

2.1 Study site and data collection

RI PT

93

The WSP is located at Ucubamba at altitude of 2,400 m a.s.l. The largest WWTP in Ecuador is

101

operated by ETAPA to treat the domestic effluent of around 400,000 inhabitants in Cuenca. The dry

102

season is between June and December and the rainy season is between January and May with

103

temperatures between 7-20 oC and 12-25 oC, respectively. This pond system is divided into two

104

identical flow lines, containing an aerated pond (AP), a facultative pond (FP) and a maturation pond

105

(MP) (Figure 1). More detailed information of Ucubamba WSP system can be found in Ho et al.

106

(2018). The data used for model calibration and validation were measured at the inlet and outlet of

107

FP1 during 2013. Following the cross validation of Snee (1977) to avoid overfitting problem of the

108

model, the data of the first half of 2013 were used for calibration while the predictive performance of

109

calibrated model were compared to the rest of the data. The area of the FP1 is 13 ha with the depth of

110

two meters, respectively, leading to the theoretical hydraulic retention time (HRT) of five days. The

111

inlet/outlet of this pond consists of a submerged pipe of 0.9 m diameter lying at the bottom of the 1.8

112

m deep pond and an overflow structure is a sharp weir of 10 m long (Alvarado et al. 2012). The

113

samples were collected during the daytime by grab sampling method three times per two weeks at the

114

input and output of the FP1, and then, were analyzed in the laboratory following the American Public

115

Health Association methods (APHA 2005). The variables were measured, including for BOD5 (mg

116

O2.L-1), COD (mg O2.L-1), ammonium (mg NH4+-N.L-1), nitrite (mg NO2--N.L-1), nitrate (mg NO3--

117

N.L-1), TKN (mg N.L-1), TP (mg P.L-1), pH (-), chlorophyll a, total solid (TS) (mg.L-1), volatile

118

suspended solid (VSS) (mg.L-1), water temperature (oC), and dissolved oxygen (DO) (mg O2.L-1).

119

Additionally, meteorological data, including air temperature (°C), solar radiation (W.m-2) and wind

AC C

EP

TE D

M AN U

SC

100

5

ACCEPTED MANUSCRIPT speed (m.s-1) were obtained from the Meteorological Station of CELEC Hidropaute located at around

121

600 m from the WSP.

122 123

Figure 1. Map of Ucubamba waste stabilization pond (WSP) in Cuenca, Ecuador.

124

2.2 Model description and implementation

125

2.2.1 System boundaries

126

A vertical 1-D model of the FP1 in Ucubamba WSP was developed by the integration of hydraulic and

127

biogeochemical submodels within the lake module of the software AQUASIM 2.1 (Reichert 1994).

128

The software is designed to simultaneously simulate physical and biogeochemical processes in natural

129

and mechanical aquatic systems. Since the pond area is relative small compared to natural lakes, the

130

assumption of homogenous horizontal gradients for the 1-D approach is satisfied (Goudsmit et al.

131

2002). As such, wind-driven horizontal flow is assumed to be negligible, hence, wind effect is

132

simulated only as the shear stress on water surface. Moreover, the bottom sediment layer is simulated

133

in a sediment submodel with an assumption of its homogeneity, meaning that no diffusive exchanges

134

among neighboring sediment layers can occur.

135

Regarding biogeochemical submodel, the growth of algae is divided into two subprocesses with

136

different nitrogen sources, NH4+ (preferred) or NO3-. Furthermore, the temporal variation of algal

137

photosynthesis and respiration activity as a result of the day-night cycle are included in the model by

138

setting the fluctuation of irradiance as a function of time. Even though light attenuation ( z) is a

139

function of algae biomass, suspended and dissolved solids, which can fluctuate within the pond

140

system, for practical purposes, it is assumed to consider

AC C

EP

TE D

M AN U

SC

RI PT

120

6

z

values as constant (Heaven et al. 2005).

ACCEPTED MANUSCRIPT The role of cyanobacteria and algal consumers is assumed to be negligible, hence being neglected. To

142

avoid complexity, only chemical equilibrium between NH4+ and NH3 is taken into account.

143

The integration of the hydraulic and biogeochemical submodels consists of a set of different

144

differential equations which are integrated numerically by DASSL algorithm follows Newton-

145

Rapson’s methods where the derivative is approximated via backward differentiation formula with

146

boundary conditions can be found in software AQUASIM 2.1 (Reichert 1994). The initial and influent

147

conditions of the model variables are supplied in Supplementary Material S1. The 1-D model was

148

implemented with the time resolution of 30 minutes, the vertical resolution of 20 cm, and the accuracy

149

was set at 0.1 %.

150

2.2.2 Hydraulic processes

151

A hydraulic submodel is built to simulate the hydrodynamic flow of the FP1, based on a combination

152

of conventional advective-diffusive model for the water column with one uniform sediment layer and

153

buoyancy-extended k-epsilon turbulence model (Elliott et al. 1999, Goudsmit et al. 2002, Rodi and

154

Research 1984, Ulrich et al. 1995). More specifically, dissolved and particulate substances are

155

transported through the water column by vertical mixing and advection. The movements of particulate

156

substances are driven by sedimentation. The bottom sediment layer is simulated in a sediment

157

submodel with an assumption of its homogeneity. The exchange of substances between water column

158

and pore water of the sediment layer via resuspension and diffusion are taken into account. The

159

surface shear is applied to count the effect of wind stress driving the horizontal motion, which is

160

estimated by following the study of Amorocho and Devries (1980).

161

In the water column, the following differential equations are solved for each dissolved and particulate

162

state variable:

165

M AN U

TE D

EP

AC C =

163 164

SC

RI PT

141

%&

=

− %&

−

+ '( # +

+ %&

−

+ ' −

⁄

−

)!* '( # +

!,

#

Eq. (1)

)+*! '!,(

Eq. (2)

where Si is the concentration of a dissolved substance i in the water column (g.m-3), Xj is the 7

ACCEPTED MANUSCRIPT 166

concentration of a particulate compound j in the water column (g.m-3), t is the time (d), z is the vertical

167

coordinate pointing downwards (z= 0 at the pond surface) (m), A is the cross-sectional area of the FP

168

(m2), Kz is the vertical mixing coefficient (m2.d-1),

169

rates in the water column (g.m-3.d-1), Q is the hydraulic flow rate (m-3.d-1), Sin and Xin are the influent

170

concentration of dissolved and particulate substances (g.m-3), hsed is the thickness of the sediment layer

171

(m), vsed is the sedimentation velocity (m.d-1), q is the volumetric lateral water flow (m2.d-1), D is the

172

coefficient of molecular diffusion (m2.d-1), θ is the porosity of the sediment layer (-), Ss,i is the

173

concentration of dissolved substance i in the sediment layer (g.m-3), Xs,j is the concentration of

174

particulate substance j in the sediment layer (g.m-3), vres is the resuspension velocity of the sediment

175

(m.d-1).

176

In the right hand side of the above equations, the first terms represent the vertical mixing because of

177

eddy diffusion, the second terms reflect the advective flow as a result of lateral influent. Especially

178

noteworthy is that the third terms represent the net effect of biogeochemical transformation processes,

179

which shows how two submodels are integrated. These transformation processes are described in more

180

detail in the biogeochemical submodel. The fourth terms are related to the pond influent. The last term

181

of Eq. (1) refers to the exchange of the dissolved substances between the water column and pore water

182

of the sediment layer. The last two terms of Eq. (2) represent the sedimentation and resuspension of

183

particulate substances between the water column and pore water of the sediment layer (Lopes et al.

184

2011, Omlin et al. 2001b). These two equations are coupled with two following equations which

185

represents the processes in the sediment layer.

187

%&

TE D

M AN U

SC

RI PT

are the biogeochemical transformation

EP

AC C

186

and

,

=

,*-,

% ,&

=

,*-,% ,& +

,

+

ℎ!* ℎ!*

,

% ,&

Eq. (3) Eq. (4)

188

In these two equations, ,*-,

189

substance j into the sediment layer (g.d-1), and

190

rates in the sediment layer (g.m-3.d-1). The first term of Eq. (3) describes the changes in mass of the

191

dissolved substances in the sediment layer due to diffusive exchange with the neighboring sediment

,

and ,*-,% ,& are the flux of the dissolved substance i and particulate ,

8

and

% &

are the biogeochemical transformation

ACCEPTED MANUSCRIPT layers if there are more than one sediment layer. The first term of Eq. (4) refers the advection of the

193

particles caused by sedimentation. The last term in both equations represents the transformation

194

processes of the substances through the sediment layer.

M AN U

195 196 197 198 199

SC

RI PT

192

Figure 2. Overview of the hydraulic processes in facultative ponds. The hydraulic submodel is based on a combination of conventional advective-diffusive model for the water column with one uniform sediment layer and buoyancy-extended k-epsilon turbulence model. The biogeochemical transformation processes were described in more detail in the biogeochemical submodel.

2.2.3 Biogeochemical processes

201

A biogeochemical submodel describes the transformations of dissolved and particulate substances

202

inside the FP1 which are specified via 30 kinetic rate equations. These mathematical equations

203

represent aerobic, anoxic, and anaerobic reactions carried out by the metabolisms and interactions of

204

algae and bacteria in the FP. The details of stoichiometric matrix, kinetic rate expressions, and their

205

values are shown in Supplementary Materials S2-S4. The formulation of the biogeochemical

206

transformation rates can be found in the Supplementary Material S5. Particularly, the aerobic and

207

anoxic processes are built following the Activated Sludge Model 2d (Henze et al. 2000), the anaerobic

208

processes are simplified from the Constructed Wetland Model No.1 (Langergraber et al. 2009) and the

209

algal metabolism is based on the River Water Quality Model 1 (Reichert et al. 2001). The Arrhenius

210

equation is included in all reaction rates to describe the temperature dependence. The interactions

211

between all involved microbial groups are summarized in Table 1.

AC C

EP

TE D

200

9

ACCEPTED MANUSCRIPT 212

Table 1. Overview of microbial competition for substrates (S) and dependencies through formed products (P) in

213

the FP.

Anaerobic populations Hydrolysis Fermenting bacteria (FB) Acetotrophic methanogens (AM) Hydrogenotrophic methanogens (HM) Anoxic and aerobic populations Anoxic heterotrophs (HB) Aerobic heterotrophs (HB) Autotrophs (AB) Algae growth on NH4+ Algae growth on NO3Phosphate accumulating organisms (PAO)

Particulate COD

Soluble COD

S

P S

Acetate

NH4+

NO3-

PO43-

O2

N2

P S

CH4

P

S S

RI PT

P

S S

H2

S S S

S S P P

P S

P

SC

S

S

P P

S

In the model, gas exchange at the pond surface includes two processes, i.e. oxygen reaeration and

215

ammonia volatilization. The oxygen reaeration rate was a flux proportional to pond depth, the

216

difference between saturation concentration and the available oxygen level, and interfacial transfer

217

coefficient (MorenoGrau et al. 1996). Concerning ammonia volatilization, in the chemical equilibrium

218

between NH4+ and NH3, the concentration of NH3 is dependent on pH, temperature, and total ammonia

219

concentration (Emerson et al. 1975). An equilibrium-based mass transfer equation in Valero and Mara

220

(2007) is applied to calculate the theoretical ammonia volatilization rate. Precipitation and

221

redissolution of phosphate PO43- are accounted, given the assumption that precipitation and

222

redissolution are reverse processes (Henze et al. 2000). Additionally, light attenuation is included to

223

describe the exponential decrease of light intensity with depth by Beer’s Law.

AC C

EP

TE D

M AN U

214

224 10

ACCEPTED MANUSCRIPT 225 226

Figure 3. Overview of biogeochemical processes in the facultative pond 1 of the Ucubamba waste stabilization

227

2.3 Parameter estimation

228

2.3.1 Local sensitivity analysis

229

Sensitivity analysis (SA) evaluates the degree to which model inputs affect model output, from that

230

the universality and robustness of these parameters can be further investigated. The following

231

technique of SA was proposed in Brun et al. (2001). Firstly, the model is defined as being described

232

above. Subsequently, the prior uncertainty of model parameters and inputs (θj) is estimated, based on

233

the literature. Scale factors of state variables, which are used to make the results for the various model

234

endpoints comparable, are calculated based on their mean concentration in the system (Brun et al.

235

2002). After prior-analysis steps, Gaussian error propagation of the uncertainty ranges is applied to

236

compute the sensitivity function for each state variable yi against the changes of model parameter θj as

237

follows. /( =

238

M AN U

SC

RI PT

pond system.

∆ & 2 !1 &

Eq. 5

where, i: 1- 7, representing the model output variables, j: 1-137, representing the model parameters,

240

∆θj represents the uncertainty ranges, which were estimated following Brun et al. (2001), Omlin et al.

241

(2001a), and Benedetti et al. (2012) ; sci is the scale factor of the state variable yi. The values of ∆θj

242

and sci are listed in Supplementary Material S6. The importance ranking of parameters was then

243

determined based on the sensitivity measure, 8(

EP

AC C

244

TE D

239

9! +

9! +

9! +

8(

, in Eq. (6)

;

= : ∑ =; / (

Eq. 6

245

A high 8(

246

while a sensitivity measure of zero indicates no effects on model output.

247

2.3.2 Identifiability analysis

248

After the calculation of the important ranking, 30 parameters with the highest values, suggesting that

249

they are the most influential parameters on the model output, was chosen for identifiability analysis.

means that the value of parameter θj has an important influence on simulation results

11

ACCEPTED MANUSCRIPT 250

Two different measures, collinearity index ? and determinant measure @, were calculated for every

251

parameter subset K of these 30 parameters to assess their identifiability. The collinearity index

252

represents the compensability of the parameter subset K which can be calculated as shown in Eq. (7). 9

;

‖C‖DE ‖!̃E GE H⋯H!̃J GJ ‖

=

9

;

‖C‖DE ‖ K G‖

=

:9

;

O K P# LMN KJ J

254

where KA being an 7 ×

255

sensitivity function /̃ A = / A ⁄‖/ A ‖, 7 is the total number of model output variables,

256

U = U; , … , UA

W

Eq. 7

RI PT

?A =

253

submatrix of the normalized matrix KA = {/̃ ( } with normalized value of the : 1-30;

is a vector of coefficient of the length K with the constraint ‖U‖ = 1; YZ is the

eigenvalue of N KAW KA P. ?A quantifies the degree of approximate linear dependence of the sensitivity

258

functions of the parameters. A value of ? over 10 indicates a poor identifiability of parameter subset

259

(Mieleitner and Reichert 2008).

260

The second criterion, determinant measure @, is defined as shown in Eq. (8).

M AN U

SC

257

;⁄ @A = [\] KAW KA #

261

A

Eq. 8

Determinant measure @ is proposed by Weijers and Vanrolleghem (1997) which combines the

263

information provided by 8(

264

value of 8(

265

2.3.3. Model calibration

266

After SA and IA, a subset of identifiable and influential parameters is chosen for calibration with the

267

data collected for the first six months of 2013. Generally, to objectify calibration process, a function

268

representing the agreement between model and data is defined to demonstrate the desire to fit the

269

model to the data (Dochain and Vanrolleghem 2001). In this case, the weighted sum of squares (WSS)

270

of the residuals is minimized from which selected parameters are calibrated using the simplex

271

algorithm (Reichert 1994). To calculate the WSS, each residual is divided by a scale factor (sci) of the

272

corresponding variable yi, from that the residuals become non-dimensional as sci and the model output

273

have a same dimension (Mieleitner and Reichert 2008).

TE D

262

9! +

, hence, good ‘‘conditional identifiability’’ of parameter subset (Brun et al. 2002).

AC C

EP

9! +

and ? in which a high value of @ indicates a low value of ? and a high

12

ACCEPTED MANUSCRIPT 2a b , c2 !1

WSS _( # = ∑ =; `

274

&#

d

Eq. 9

In this equation, f9*g!, is the i-th measurement, f _( # is the calculated value of the model output

276

corresponding to the i-th measurement and evaluated at the time and location of this measurement, and

277

n is the number of data points.

278

2.4 Uncertainty analysis

279

Like natural lakes, WSP systems are inherently subject to large variations, hence, uncertainty analysis

280

(UA) is considered an essential role during their modeling process (Ho et al. 2017). Since SA can help

281

to recognize the uncertainty source from model inputs which allocates to the uncertainty of model

282

output, it is recommended that UA and SA should be implemented in tandem (Saltelli et al. 2008). In

283

this study, UA is conducted in two cases. In the first case, all of the prior uncertainties listed in

284

Supplementary Material S6 are included while the calibrated value of the identifiable parameters is

285

fixed with no uncertainty in the second case. These calculations allow assessing the contribution of

286

these parameters to the uncertainty of model output. Monte Carlo simulation is applied using the

287

UNCSIM package (Reichert 2005) to propagate the model uncertainty into the predictions (Mckay et

288

al. 1979). 500 sets of samples are generated using Latin Hypercube Sampling (LHS) technique, which

289

provides a sufficient coverage of parameter space.

290

3. Results and discussion

291

3.1 Sensitivity analysis

292

Sensitivity analysis (SA) is designed as a tool to identify the most influential model parameters for the

293

variability of the state variables. For the overparameterized models, SA is considered very useful as

294

model output are often strongly influenced by few key inputs (Saltelli 2002). In this case study, we

295

investigate particularly the degree to which model inputs affecting two groups of model output, i.e.

296

organic matter and nutrient removals.

297

3.1.1 The most influential parameters for organic matter removal

298

Organic matter (OM) content in wastewater can be measured via COD which includes slowly

AC C

EP

TE D

M AN U

SC

RI PT

275

13

ACCEPTED MANUSCRIPT biodegradable particulate COD (XS), fermentable and readily biodegradable soluble COD (SF),

300

fermentation products as acetate (SA), inert soluble (SI) and particulate (XI) COD, and COD from

301

bacterial biomass. However, not all of these components are of equal importance (Henze et al. 2000).

302

Particularly, COD fraction from microorganisms is mostly neglected due to wash-out phenomena in

303

treatment systems. Additionally, the inert COD is not in our interest because of its marginal variation.

304

To identify the most influential model parameters for the first three COD fractions, the sensitivity

305

functions (/ ( ) for each of the three variables against the changes of model parameter θj were

306

calculated. Positive values of / ( indicate the proportional relationship between the model parameters

307

and the concentration of state variable and vice versa. As shown in Figure 4, the most significant

308

parameter influencing the OM content is the coefficient of vertical turbulent diffusion (Kz) which is

309

responsible for around 60%, 40%, and 20% of the total variance of XS, SF, and SA, respectively.

310

Interestingly, while the increase of Kz has a positive effect on the accumulation of XS, meaning a

311

negative impact on its removal, higher Kz leads to higher removal efficiency of SF and SA. This result

312

suggests that high vertical mixing leading to non-stratified water column is detrimental to hydrolysis, a

313

conversion process of XS to SF conducted by heterotrophs (XH) and fermenting bacteria (XFB). Indeed, a

314

similar oxygen level in water column can inhibit the growth of anaerobic bacteria at bottom layers,

315

like XFB and anoxic heterotrophs (XH,NO) in this case, while providing an insufficiently amount of

316

oxygen for aerobic heterotrophs (XH,O) at the surface layers. The removal of XS via hydrolysis creating

317

SF and then SA via the growth of XFB in the wastewater explain the abovementioned contradictory

318

effects of the Kz. Noticeably, the maximum growth rate of heterotrophs (µ H) is responsible for 20% of

319

the total variance of SF and SA. This parameter is highly dependent on temperature as the optimal

320

growth rate of XH in WSPs was reported at 35 oC, being similar to the recommended value for

321

nitrifying bacteria (Mayo and Noike 1996). Indeed, temperature appears to be a crucial factor for the

322

efficiency of OM removal with high contribution of temperature coefficient (θTw) in all graphs.

AC C

EP

TE D

M AN U

SC

RI PT

299

14

323 324 325 326 327

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 4. Ten most significant parameters influencing the variability of slowly biodegradable particulate COD (XS), fermentable and readily biodegradable soluble COD (SF), and fermentation products as acetate (SA). Besides inert fractions, these variables represent nearly all of organic matter content in municipal wastewater. The values in the radar graphs are the results of sensitivity function (/ ( ) for each state variable against the changes of model parameters.

329

3.1.2 The most influential parameters for nutrient removal

330

Nutrient content in municipal wastewater mainly includes two fractions of nitrogen (N) and

331

phosphorus (P). Particularly, the total N concentration contains particulate and soluble Kjeldahl

332

nitrogen (TKN), and nitrate- and nitrite-nitrogen (SNO) (Henze et al. 2000). Particulate Kjeldahl

333

nitrogen, as the sum of N bound to all organic particulate fractions (XI, XS, XH, XPAO, and XA), is not

334

investigated due to its marginal values compared to the value of soluble Kjeldahl nitrogen which is

335

dominated by ammonium-nitrogen (SNH). For the same reason, particulate phosphorus is not in the

336

interest of this research. As such, the most influential model parameters for nutrient removal

337

represented by the variance of SNH, SNO, and SPO4 are shown in Figure 5. Regarding the first state

338

variable, it appears that microbial assimilation and ammonia volatilization are two main processes of

339

SNH removal in FPs. The first process is demonstrated via the role of the metabolisms of algae,

340

autotrophic (XA) and methanogenic microorganisms (XAM and XHM). Specifically, the decay (bALG) and

AC C

EP

TE D

328

15

ACCEPTED MANUSCRIPT respiratory rates (h +*! ij of algae are responsible for above 20% of the total ammonium removal which

342

is equal to the total contribution of the bacterial metabolisms on ammonium uptake. More importantly,

343

pH contributes to more than 20% of the total variance of SNH, demonstrating ammonia volatilization as

344

one of the most important processes for ammonium removal. In fact, ammonia volatilization can

345

account for 75-98% of total N removal with pH range from 7 to 9 (Pano and Middlebrooks 1982,

346

Pearson et al. 1996, Reed 1985). These results indicated both advantages and disadvantages of the

347

unique characteristic of high-altitude climate in this case. On the one hand, strong solar radiation

348

boosted the algae photosynthesis, leading to more nitrogen algal uptake and high value of pH, which

349

can increase the volatilization rate. On the other hand, the abundance of algae can release the

350

ammonium back to the water column during their decay process, which is indicated via the positive

351

correlation between the algal decay rate and SNH concentration.

352

Due to the limited concentration of SNO from the influent of the FP1, its main input is from nitrification

353

process by autotrophic bacteria. As such, the metabolism of this microorganism, represented as

354

maximum growth rate (µ A), yield coefficient (YA), and phosphorous half saturation coefficient (

355

the radar chart, contributes above 50% of the total variance of SNO. The absence of parameters related

356

to heterotrophic bacteria indicates the insignificant role of denitrification process in the removal of

357

SNO, which can be explained by the aerobic conditions in the pond. On the other hand, the Kz is

358

responsible for around 15%, suggesting the stratification in water column is crucial for the presence of

359

autotrophs as they require oxygen generated from algal photosynthesis.

360

The following parameters were found to have a significant impact on the P removal (listed in the

361

decreasing order of importance): fraction of phosphorus in SF (iPSF), yield coefficient of fermenting

362

bacteria (YFB), decay rate of algae (bALG), and fraction of P in bacteria (iPSB). This list suggests the

363

negligible role of chemical precipitation process in P removal in Ucubamba WSP systems. In fact, the

364

possibility of its occurrence is constrained by marginal concentration of cations (Al3+ and Fe3+),

365

around only 4 mg.L-1 while the sufficient amount is suggested larger than 50-100 mg.L-1 (Diaz et al.

366

1994). As a result, the main process of P removal is the uptake by microorganisms which, however,

367

has a marginal effect on the P removal, leading to its similar concentration in the effluent.

k

) in

AC C

EP

TE D

M AN U

SC

RI PT

341

16

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

368 369 370 371 372

Figure 5. Ten most significant parameters influencing the variability of ammonium-nitrogen (SNH), nitrate- and

373

3.1.3 The most important parameters for model output

374

After calculating / ( , further step of SA is proceeded with the computation of importance rankings

375

(8 9! + ) of model parameters. Table 2 shows the 8 9!

376

on the model output. The ranking suggests that the physical parameters, i.e. eddy vertical diffusivity

377

(Kz), water temperature (Tw), and temperature coefficient (θTw), initiate the most variability of the

378

model output. Moreover, the metabolisms of autotrophic and heterotrophic bacteria, indicated via

379

maximum specific growth and yield coefficient (µ A, YA, YH), also appear as highly influential

380

parameters. In fact, autotrophs prove their significance by being responsible for 50% of the total

381

variance of SNO and µ H is one of the most influential factors on OM removal. More importantly, high

382

8 9!

383

output towards the efficiency of the photosynthetic process. Indeed, the aeration of FPs depends

384

heavily on algal activity where more than 80% of dissolved oxygen is generated in WSP system

nitrite-nitrogen (SNO), and inorganic phosphorus (SPO4). Besides particulate fractions, these variables represent nearly all of nutrient content in municipal wastewater. The values in the radar graphs are the results of sensitivity

+

of 30 parameters which have highest influence

AC C

EP

TE D

function (/ ( ) for each state variable against the changes of model parameters.

+

of algal-related parameters (pH, µ ALG, kz, h +*! ij , fp1, bALG) reflect high sensitivity of the model

17

ACCEPTED MANUSCRIPT 385

(Shilton 2005). This is certainly the case in this high-altitude WSP where strong solar radiation, up to

386

1500 W.m-2 during the peak of light intensity, can promote the overgrowth of algae, generating

387

extremely high oxygen level up to 39 mg O2.L-1 (Ho et al. 2018). This prevailing of aerobic condition

388

reduce the influence of anaerobic processes on the performance of the pond system, which reflects in

389

low 8 9!

390

Table 2. Importance rankings (8 9! + ) of the top 30 parameters

8 9! + 4.025 2.267 0.444 0.429 0.428 0.292 0.252 0.196 0.193 0.191

Rank 11 12 13 14 15 16 17 18 19 20

Parameters kz h +*! ij

YPHA YFB

8 9! + 0.114 0.066 0.059 0.038 0.037 0.036 0.032 0.031 0.021 0.018

Rank 21 22 23 24 25 26 27 28 29 30

Parameters YPO4 YAM bH bHM

SC

Parameters Kz pH Tw YPAO θTw µA YA µH µ ALG YH

RI PT

of anaerobic-related parameters (YFB, YHM, Kh, YAM, bHM).

l

fp1 iPSF bALG YHM Kh

M AN U

Rank 1 2 3 4 5 6 7 8 9 10

+

n m n l

Ta iPXB nP iNSF

8 9! + 0.018 0.016 0.013 0.012 0.011 0.010 0.006 0.006 0.006 0.005

Physical-related parameters: Kz: eddy vertical diffusivity; Tw: water temperature; θTw: temperature coefficient; Ta: air temperature. Anaerobic-related parameters: YFB: Yield of FB; YHM: Yield of HM; Kh: Hydrolysis rate; YAM: Yield of AM; bHM: decay rate of HM. Heterotrophs-related parameters: µ H: Maximum growth of heterotrophs; YH: Yield of heterotrophs; bH: decay rate of heterotrophs; mn : Fermentable substrate half saturation of heterotrophs; ln : Oxygen half saturation of heterotrophs. Autotrophs-related parameters: µ A: Maximum growth of autotrophs; YA: Yield of autotrophs; +*! l : Oxygen half saturation of autotrophs. Algae-related parameters: pH; µ ALG: Maximum growth of algae; h ij : Respiration rate of algae; fp1: Fraction of XI formed during the decay of algae; bALG: Decay rate of algae. PAO-related parameters: YPAO: Yield of PAO; YPHA: PHA requirement for PP storage; iPSF: fraction of P in XS; YPO4: PP requirement (SPO4 release) for PHA storage; iPXB: Fraction of P in bacteria; nP: Correction factor for µ PAO under anoxic conditions

400

3.2 Identifiability analysis

401

To avoid an inevitable problem of poor identifiability of the complex models, the compensability of

402

model parameters is evaluated via collinearity index ? and determinant measure @. Figure 6 illustrates

403

an overview of possible values of the two measures for all the possible subsets of the 30 parameters

404

listed in Table 2. i.e. 2.65e+32 combinations. In contrast to the dramatic drop of ρ, the graph

405

demonstrates the significant increase of the maximum values of ? with the increasing size of the

406

parameter subsets. Interestingly, the highest ? value of the subset of two parameters already increases

407

to 28.01, suggesting that a set of two parameters can have a serious problem associated with parameter

408

correlation as the threshold value is 10 (Omlin et al. 2001a). This wide range of ? also ascertains the

409

necessity of IA in model evaluation to obtain the uniqueness of parameter estimates and to avoid the

410

poor convergence in nonlinear parameter estimation (Vanrolleghem et al. 1995). Especially

AC C

EP

TE D

391 392 393 394 395 396 397 398 399

18

ACCEPTED MANUSCRIPT noteworthy is that when the size of subset is above 8, the value range of @ diminishes considerably,

412

which can make the differentiation of the results among different subsets difficult.

M AN U

413 414 415 416

SC

RI PT

411

Figure 6. Overview of the ranges of collinearity index (?) and determinant measure (@) for all possible parameter subsets of all size (1-30 parameters). The blue lines indicate the area when serious problems of identifiability start, i.e. ? between 10 and 15 (Reichert and Vanrolleghem 2001).

To avoid assessing the collinearity of all combinations, the parameters were divided into five different

418

functional subsets, based on the mechanisms of the FP. The reason of this division is that the

419

parameters in a same functional subset have higher chance to represent a similar property of the

420

system, hence, these parameters are more likely to have high correlation. By identifying them and then

421

removing them from the final combination, the simultaneous estimation of this combination is more

422

likely to be successful. As such, five functional subsets are defined, i.e. physical, anaerobic, aerobic

423

and anoxic, algal, and PAO-related processes, to investigate their correlations via ? and @ (see Table

424

3). The details of the parameters belonging to the five subsets can be found in the Supplementary

425

Material S7. The first group evidences high correlation between water temperature (Tw) and

426

temperature coefficient (θTw) as adding θTw into the subset increases the ? by 12. On the other hand, the

427

parameters related to anaerobic process and algal activity shows relatively low correlations. While

428

three anaerobic parameters (YFB, YHM, Kh) describing different anaerobic processes explain their weak

429

interdependence, low correlations among algal-related parameters are surprising outcomes.

430

Theoretically, there is a strong relation between the algal activity and pH since CO2 consumption via

AC C

EP

TE D

417

19

ACCEPTED MANUSCRIPT algal photosynthesis causes the increase of pH during the day while its respiration generates CO2,

432

causing the drop of pH and O2 level during the night or light limited conditions. Since the model

433

excludes the chemical equilibria between CO2 and HCO3-, between HCO3-and CO32-, between H2O

434

and H+ and OH- due to its current complexity, this carbonate bicarbonate buffering system in FPs is

435

not able to be included. Regarding aerobic and anoxic process, the maximum growth rate (µ A) and

436

yield coefficient (YA) of autotrophs have a high correlation, above 6, while the connection between

437

autotroph- and heterotroph-related parameters is insignificant.

438

After identifying the set of parameters having strong interdependency, all possible combinations of the

439

parameters from the five functional group excluding these sets were calculated. However, since the

440

obtained lowest ? of the combinations with more than 13 parameters was already 13.42, implying the

441

identifiability problem, the ? and ρ are calculated to the possible combinations of size 12 and 13.

442

These combinations can include four physical parameters, three anaerobic parameters, five algal-

443

related parameters, five parameters link to autotrophs and heterotrophs, and three PAO-related

444

parameters. In contrast to similar values of @, ? changes considerably among different combined

445

subsets, from 8 to 29, suggesting its potential of being a more relevant criterion in assessing the

446

identifiability of the model parameters. As such, after assessing both indices of the combinations, the

447

subset in bold in Table 3 is selected based on its lowest value of ? and relatively high value of ρ

448

compared to other subsets. The highest values of ? are observed in the subsets containing two

449

parameters, i.e. µ ALG and light extinction coefficient (kz), suggesting the identifiability problem from

450

the parameterization of light dependence of algal growth which is line with the conclusions of Omlin

451

et al. (2001a) on a natural lake system. In fact, the intrinsic correlation between algal biomass and

452

light extinction coefficient has been studied and modelled using a (non)linear approximation in

453

photobioreactors (Fernandez et al. 1997, Ogbonna et al. 1995, Privoznik et al. 1978, Yun and Park

454

2001). This suggests that there is a need for further investigation with more available data, that can

455

indicate the spatiotemporal variability of algal biomass in the WSPs. From that, a (non)linear equation

456

between this parameter and kz could be added to improve the model structure, hence, its identifiability.

AC C

EP

TE D

M AN U

SC

RI PT

431

20

ACCEPTED MANUSCRIPT 457

Table 3. Collinearity index (?) and determinant measure (@) of selected parameters subsets. Based on their

458

function, the parameters are categorized into five process groups. The subset in bold is selected for calibration

459

based on the values of ? and @. The description of the parameters can be found in the footnote of Table 2.

Aerobic & Anoxic

PAO

Combination

Parameters Kz, Tw, kz Kz, Tw, θTw Kz, Tw, θTw, kz YFB, YHM YFB, Kh YFB, YHM, Kh pH, µ ALG, fp1 pH, µ ALG, bALG pH, µ ALG, h +*! ij pH, µ ALG, fp1, bALG pH, µ ALG, h +*! ij , fp1 µ A, YA, µ H µ A, YA, YH µ A, YA, l µ A, YA, µ H, YH µ A, YA, µ H, l YPAO, YPHA YPAO, iPSF YPAO, YPHA, iPSF Kz, Tw, YFB, pH, µALG, µA, Kh, YPAO, YPHA, YHM, YA, YH Kz, Tw, YFB, pH, µ ALG, µ A, Kh, YPAO, YPHA, bALG, YA, YH Kz, Tw, YFB, pH, µ ALG, µ A, Kh, YPAO, YPHA, µ H, YH, iPSF Kz, Tw, YFB, pH, µ ALG, µ A, Kh, YPAO, YPHA, fp1, µ H, YA Kz, Tw, YFB, pH, µ ALG, µ A, Kh, YPAO, YPHA, bALG, µ H, kz Kz, Tw, YFB, pH, µ ALG, µ A, Kh, YPAO, YPHA, fp1, µ H, YH, kz

RI PT

Algal activity

@ 14.74 11.21 7.91 0.72 0.67 0.62 6.31 5.55 5.53 3.49 3.15 3.69 3.66 3.29 3.86 1.80 2.42 3.02 1.69 3.47 3.37 3.45 3.19 2.51 2.75

SC

Anaerobic

? 1.24 12.12 12.71 1.13 1.01 1.13 1.20 1.72 3.94 1.84 6.32 6.78 6.80 7.32 6.80 7.34 3.96 1.01 3.97 8.16 8.31 13.43 18.73 29.04 29.06

M AN U

Physical

Set size 3 3 4 2 2 3 3 3 3 4 4 3 3 3 4 4 2 2 3 12 12 12 12 12 13

TE D

Group

3.3 Model calibration

461

After the establishment of the subset of 12 selected parameters, an automatic calibration was

462

performed on the basis of the measurements until the convergence in simplex mode is achieved in

463

AQUASIM (Reichert 1994). The original and calibrated results are represented in Table 4. The first

464

conclusion can be drawn is that, except for YA, stoichiometric parameters maintain comparable values

465

after the calibration, inferring their less influential role on the model output compared to the kinetic

466

parameters. Conversely, the maximum growth of algae (µ ALG) increase considerably by 49.1%

467

indicating more intensive algal activity in this pond. In fact, Ucubamba WSP locating at 2,400 m a.s.l

468

receives immense light intensity, hence, accelerating the photosynthesis rate and the growth of algae

469

(Ho et al. 2018). Higher calibrated pH of 9.09 can also indicate this unique characteristic of high

470

altitude pond where volatilization process was one of the main nitrogen removal processes within the

AC C

EP

460

21

ACCEPTED MANUSCRIPT pond. However, note that since the carbonate bicarbonate buffering system in the FP1 was not

472

included in the model, the drop of pH during the night was not able to be simulated. This might lead to

473

the overestimation of the rate of volatilization process during the night. Two physical parameters (Kz

474

and Tw) experience a dramatic increase after the calibration. As the field data were collected during the

475

day, Tw would be higher than the daily average, 18.6 oC while extreme high-altitude climates, i.e.

476

strong solar radiation, great temperature variation and strong wind, can be a reason for the higher

477

calibrated value of Kz. Compared to the ordinary range of vertical eddy diffusivity in shallow and

478

eutrophic lakes, from 10-7 to 10-4 m2.s-1 (Macintyre 1993), the calibration presented a relatively high

479

value of Kz, 4.86 ×10-5 m2.s-1. Since Kz characterizes the intensity of vertical mixing in the FP1, its

480

high value indicates the homogenous conditions of the water column in the FP1 due to abundant

481

aerobic condition and strong mixing. The aerobic prevailing condition in Ucubamba pond system

482

during the daylight is also implied via a substantial drop in value of hydrolysis rate (Kh) by almost

483

50%. Furthermore, the reduction of nitrification, reflecting by the decrease of µ A and YA by 17.8% and

484

49.1%, respectively, supports the findings of Craggs et al. (2000) who concluded that the lack of

485

biofilm attachment surfaces is a major constraint for its efficiency in WSPs. As such, the most

486

influential parameters for the N removal are mainly related to microbial assimilation and ammonia

487

volatilization.

488

Table 4. Original and calibrated values of the selected parameters.

EP Unit m2.d-1 d-1 d-1 d-1 o C g COD.g-1 N g COD.g-1 COD g COD.g-1 COD g COD.g-1 COD g COD.g-1 COD g COD.g-1 COD

AC C

Parameter Kz Kh µA µ ALG pH Tw YA YFB YH YHM YPAO YPHA

TE D

M AN U

SC

RI PT

471

Original value 0.05 3 2 2 7.8 18.6 0.24 0.053 0.63 0.02 0.63 0.2

Calibrated value 0.07 1.53 1.64 2.98 9.09 25.1 0.12 0.05 0.60 0.02 0.66 0.21

Change (%) 49.3 -49.0 -17.8 49.1 16.6 34.8 -49.1 -4.6 -4.8 -4.9 4.7 4.9

489

Table 5 shows better predictive performance of the calibrated model in estimating the removal

490

efficiencies of the FP via relatively lower values of two error measurements, i.e. WSS and mean 22

ACCEPTED MANUSCRIPT absolute percentage error (MAPE). By using these two error measurements, the model goodness of fit

492

with respect to mean, amplitude and phase can be accessed. In fact, except for BOD, the model

493

predicts relatively precisely the removal of OM and nutrients. Its inaccuracy in predicting BOD

494

concentration can be due to the fact that the biogeochemical submodel is built on COD balance. The

495

changes in the value of parameters related to algal photosynthesis and physical processes appear to

496

reduce considerably the error measures of the model, i.e. 42.2 % and 54.7 % of WSS of BOD and TN.

497

Table 5. Weighted sum of squares (WSS) and mean absolute percentage error (MAPE) for initial and end values

498

(after calibration). MAPE formula can be found in Supplementary Material S8. WSSend 90.1 33.2 3.1 7.4

MAPEsts (%) 65.9 27.3 14.9 10.7

MAPEutv (%) 48.9 33.4 9.0 10.2

SC

WSSini 155.9 31.2 6.9 7.0

Change (%) 42.2 -6.3 54.7 -5.6

M AN U

Parameters BOD COD TN TP

RI PT

491

Change (%) 25.8 -22.2 39.6 4.6

3.4 Uncertainty analysis

500

Figure 7 shows the uncertainty range of model predictions, i.e. OM and nutrient removals, in

501

comparison with the observed data collected during the second half of 2013. It is important to

502

recognize the substantial contribution of the selected parameters on the output uncertainty ranges

503

since, apart from the case of N removal, the light grey bands are much wider than the dark grey bands.

504

This suggests that the consideration of model calibration can reduce greatly the model uncertainty and

505

the selected parameters appear to be key inputs influencing the model output. However, these

506

relatively small uncertainty ranges with no contribution of the selected parameters are not able to

507

encompass the large variations of measured data. This high variabilities of the effluent concentrations

508

of FP1 can be associated with the large variations of climatic conditions which are very influential

509

factors on the performance of high-altitude pond systems (Ho et al. 2018). On the other hand, a wide

510

dark-grey range in the prediction of N removal infers that there is/are other parameter(s) besides the

511

selected parameters, which can largely contribute to the model uncertainty. Also noteworthy is the

512

extensive light-grey uncertainty ranges, which indicate this intrinsic property of large model

513

predictions, especially in the case of ecological systems such as ponds, lakes, and rivers (Omlin et al.

514

2001a, Reichert and Vanrolleghem 2001). Regarding the model predictive performance, besides the

AC C

EP

TE D

499

23

ACCEPTED MANUSCRIPT inaccuracy in predicting BOD via its high overestimation, the model shows relatively accurate

516

predictions for COD and nutrient removal, with minor underestimations. These underestimations can

517

be caused by the seasonal difference in sampling period of the data used in the calibration and

518

validation processes. Particularly, the measurements used for calibration were collected during the

519

first half of 2013, which was a rainy season while the validated data were measured during the dry

520

season. As such, the dilution factor from heavy rains may cause the underestimations of model

521

predictions in dry season. More importantly, the patterns of collinearity among model parameters,

522

such as temperature and its coefficients or algae growth rate and light attenuation, are likely to change

523

in different environmental conditions. Additionally, the difference regarding the characterization of

524

the influent between two seasons might also another reason of the inaccuracy of the model. From this

525

perspective, a well-designed sampling strategy using online sensor technology, that is able to cover

526

representative geographic and environmental conditions, such as night-time data or different seasonal

527

data, is needed for the model extrapolation to longer time period and to other pond treatment systems.

AC C

EP

TE D

M AN U

SC

RI PT

515

24

528 529 530 531 532 533 534

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 7. Uncertainty range of model predictions, i.e. organic matter and nutrient removal efficiencies, in comparison with the observed data collected during the second half of 2013. The dashed lines represent mean predicted results while the dot-dashed lines represent the observed data. The light grey bands represent the uncertainty range using 10th and 90th percentile values of the Monte Carlo simulations with all prior uncertainties while the dark grey bands represent the uncertainty range without the contribution of the selected parameters, which are fixed with the calibrated values in Table 4.

25

ACCEPTED MANUSCRIPT 4. Conclusions

536

• We develop an all-encompassing mechanistic model simulating not only biogeochemical processes

537

of carbon, nitrogen, and phosphorus removal but also the hydraulic and physical processes of a

538

facultative pond in combination with a systematic investigation for model calibration and

539

validation. Despite the high variability of the measurements, the model shows a good agreement

540

between the predictions and observed data, indicating via the low values of two error

541

measurements, i.e. WSS and MAPE.

RI PT

535

• The first application of a practical sensitivity analysis sheds light on the most influential parameters

543

on the pond performance. Particularly, together with the maximum growth rate of heterotrophs,

544

physical parameters, such as eddy vertical diffusivity, water temperature, and temperature

545

coefficient, initiate the most variability of organic matter. Nitrogen removal is mainly done by

546

microbial assimilation and ammonium volatilization in contrast to insignificant role of

547

denitrification process due to the abundant aerobic condition in the high-altitude pond as a result of

548

intensive algae photosynthetic activity. Phosphorus removal by phosphorus-accumulating

549

organisms and chemical precipitation appears to be trivial, leading to low removal efficiency.

TE D

M AN U

SC

542

• Two indexes are applied for identifiability analysis, i.e. collinearity index ? used by Brun et al.

551

(2001) and determinant measure @ proposed by Weijers and Vanrolleghem (1997). Based on their

552

values of all possible parameter subsets, the former appears to be a more relevant indicator for

553

assessing the parameter identifiability, while the latter emerges as an excellent tool for determining

554

the most influential parameter subsets.

AC C

EP

550

555

• A deeper understanding of parameter identifiability problems can be obtained thanks to the

556

practical identifiability analysis. Particularly, two main problems are identified. The first problem

557

comes from the collinearity of maximum growth rate of algae and light attenuation coefficient

558

indicating the strong influence of algal biomass on the light penetration in the pond system. To

559

reduce the problem of this collinearity, better understanding about the mechanisms of this effect

560

should be added to improve the model structure, hence, its identifiability. Secondly, like the pond

561

performance, the patterns of collinearity among model parameters are likely to change in different

26

ACCEPTED MANUSCRIPT environmental conditions, such as temperature and its coefficients in this case, causing serious

563

errors in model extrapolation. As such, a well-designed sampling strategy using online sensor

564

technology, that is able to cover representative geographic and environmental conditions, such as

565

night-time data or different seasonal data, is needed for the model extrapolation to more extensive

566

time period and to other pond treatment systems.

RI PT

562

• Model calibration reveals the distinctive characteristics of this high-altitude pond where more

568

intensive algal activity because of strong solar radiation generating aerobic prevailing condition.

569

Consequently, this leads to low hydrolysis rate, low importance ranking of anaerobic-related

570

parameters, and the absence of anoxic heterotrophic bacteria which are responsible for

571

denitrification process. Higher values of eddy vertical diffusivity and water temperature are

572

obtained, reflecting the impact of high-altitude climates, i.e. great variation of temperature and

573

strong wind.

M AN U

SC

567

• An uncertainty analysis suggests that a large avoidable uncertainty as a result of vast complexity of

575

the applied model can be reduced greatly via a systematic procedure for parameter estimation.

576

However, one should keep in mind the fluctuation of the model parameters when extrapolating the

577

model output, especially in this case where seasonal variation and the characterization of the

578

influent appear to be influential factors.

TE D

574

Acknowledgement

580

This research was performed in the context of the VLIR Ecuador Biodiversity Network project. This

581

project was funded by the Vlaamse Interuniversitaire Raad-Universitaire Ontwikkelingssamenwerking

582

(VLIR-UOS), which supports partnerships between universities and university colleges in Flanders

583

and the South. We are grateful to ETAPA for providing the data. Long Ho is supported by the special

584

research fund (BOF) of Ghent University.

585

Declarations of interest

586

None

AC C

EP

579

27

ACCEPTED MANUSCRIPT References

588 589 590

Alvarado, A., Sanchez, E., Durazno, G., Vesvikar, M. and Nopens, I. (2012) CFD analysis of sludge accumulation and hydraulic performance of a waste stabilization pond. Water Science and Technology 66(11), 2370-2377.

591 592

Amorocho, J. and Devries, J.J. (1980) A New Evaluation of the Wind Stress Coefficient over Water Surfaces. Journal of Geophysical Research-Oceans 85(Nc1), 433-442.

593 594

APHA (2005) Standard methods for the examination of water and wastewater. American Public Health Association (APHA): Washington, DC, USA.

595 596 597

Benedetti, L., Batstone, D.J., De Baets, B., Nopens, I. and Vanrolleghem, P.A. (2012) Uncertainty analysis of WWTP control strategies made feasible. Water Quality Research Journal of Canada 47(1), 14-29.

598 599

Beran, B. and Kargi, F. (2005) A dynamic mathematical model for wastewater stabilization ponds. Ecological Modelling 181(1), 39-57.

600 601

Brun, R., Kuhni, M., Siegrist, H., Gujer, W. and Reichert, P. (2002) Practical identifiability of ASM2d parameters - systematic selection and tuning of parameter subsets. Water Research 36(16), 4113-4127.

602 603

Brun, R., Reichert, P. and Kunsch, H.R. (2001) Practical identifiability analysis of large environmental simulation models. Water Resources Research 37(4), 1015-1030.

604 605

Colomer, F.L. and Rico, D.P. (1993) Mechanistic Model for Facultative Stabilization Ponds. Water Environment Research 65(5), 679-685.

606 607 608

Craggs, R.J., Tanner, C.C., Sukias, J.P.S. and Davies-Colley, R.J. (2000) Nitrification potential of attached biofilms in dairy farm waste stabilisation ponds. Water Science and Technology 42(10-11), 195-202.

609 610

Diaz, O.A., Reddy, K.R. and Moore, P.A. (1994) Solubility of Inorganic Phosphorus in Stream Water as Influenced by Ph and Calcium-Concentration. Water Research 28(8), 1755-1763.

611 612

Dochain, D., Gregoire, S., Pauss, A. and Schaegger, M. (2003) Dynamical modelling of a waste stabilisation pond. Bioprocess and Biosystems Engineering 26(1), 19-26.

613 614

Dochain, D. and Vanrolleghem, P.A. (2001) Dynamical Modelling & Estimation in Wastewater Treatment Processes, IWA Publishing.

615 616

Elliott, J.A., Irish, A.E., Reynolds, C.S. and Tett, P. (1999) Sensitivity analysis of PROTECH, a new approach in phytoplankton modelling. Hydrobiologia 414, 45-51.

617 618 619

Emerson, K., Russo, R.C., Lund, R.E. and Thurston, R.V. (1975) Aqueous Ammonia Equilibrium Calculations - Effect of Ph and Temperature. Journal of the Fisheries Research Board of Canada 32(12), 2379-2383.

620 621 622

EPA, U. (2011) Principles of Design and Operations of Wastewater Treatment Pond Systems for Plant Operators, Engineers, and Managers, United States Environmental Protection Agency, Office of Research and Development.

623 624 625

Fernandez, F.G.A., Camacho, F.G., Perez, J.A.S., Sevilla, J.M.F. and Grima, E.M. (1997) A model for light distribution and average solar irradiance inside outdoor tubular photobioreactors for the microalgal mass culture. Biotechnology and Bioengineering 55(5), 701-714.

AC C

EP

TE D

M AN U

SC

RI PT

587

28

ACCEPTED MANUSCRIPT Gehring, T., Silva, J.D., Kehl, O., Castilhos, A.B., Costa, R.H.R., Uhlenhut, F., Alex, J., Horn, H. and Wichern, M. (2010) Modelling waste stabilisation ponds with an extended version of ASM3. Water Science and Technology 61(3), 713-720.

629 630 631

Goudsmit, G.H., Burchard, H., Peeters, F. and Wuest, A. (2002) Application of k-epsilon turbulence models to enclosed basins: The role of internal seiches. Journal of Geophysical Research-Oceans 107(C12).

632 633

Heaven, S., Banks, C.J. and Zotova, E.A. (2005) Light attenuation parameters for waste stabilisation ponds. Water Science and Technology 51(12), 143-152.

634 635

Henze, M., Guijer, W., Mino, T. and van Loosdrecht, M. (2000) Activated sludge models ASM1, ASM2, ASM2d and ASM3, IWA Publishing, IWA Publishing.

636 637 638

Ho, L., Pham, D., Van Echelpoel, W., Muchene, L., Shkedy, Z., Alvarado, A., Espinoza-Palacios, J., Arevalo-Durazno, M., Thas, O. and Goethals, P. (2018) A Closer Look on Spatiotemporal Variations of Dissolved Oxygen in Waste Stabilization Ponds Using Mixed Models. Water 10(2), 201.

639 640

Ho, L.T., Van Echelpoel, W. and Goethals, P.L.M. (2017) Design of waste stabilization pond systems: A review. Water Research 123, 236-248.

641 642 643

Kayombo, S., Mbwette, T.S.A., Mayo, A.W., Katima, J.H.Y. and Jorgensen, S.E. (2000) Modelling diurnal variation of dissolved oxygen in waste stabilization ponds. Ecological Modelling 127(1), 2131.

644 645 646

Langergraber, G., Rousseau, D.P.L., Garcia, J. and Mena, J. (2009) CWM1: a general model to describe biokinetic processes in subsurface flow constructed wetlands. Water Science and Technology 59(9), 1687-1697.

647 648 649

Lopes, F., Viollier, E., Thiam, A., Michard, G., Abril, G., Groleau, A., Prevot, F., Carrias, J.F., Alberic, P. and Jezequel, D. (2011) Biogeochemical modelling of anaerobic vs. aerobic methane oxidation in a meromictic crater lake (Lake Pavin, France). Applied Geochemistry 26(12), 1919-1932.

650 651

Macintyre, S. (1993) Vertical Mixing in a Shallow, Eutrophic Lake - Possible Consequences for the Light Climate of Phytoplankton. Limnology and Oceanography 38(4), 798-817.

652 653

Mayo, A.W. and Noike, T. (1996) Effects of temperature and pH on the growth of heterotrophic bacteria in waste stabilization ponds. Water Research 30(2), 447-455.

654 655 656

Mckay, M.D., Beckman, R.J. and Conover, W.J. (1979) A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code. Technometrics 21(2), 239-245.

657 658

Mieleitner, J. and Reichert, P. (2008) Modelling functional groups of phytoplankton in three lakes of different trophic state. Ecological Modelling 211(3-4), 279-291.

659 660 661

MorenoGrau, S., GarciaSanchez, A., MorenoClavel, J., SerranoAniorte, J. and MorenoGrau, M.D. (1996) A mathematical model for waste water stabilization ponds with macrophytes and microphytes. Ecological Modelling 91(1-3), 77-103.

662 663

Ogbonna, J.C., Yada, H. and Tanaka, H. (1995) Kinetic-Study on Light-Limited Batch Cultivation of Photosynthetic Cells. Journal of Fermentation and Bioengineering 80(3), 259-264.

664 665

Omlin, M., Brun, R. and Reichert, P. (2001a) Biogeochemical model of Lake Zurich: sensitivity, identifiability and uncertainty analysis. Ecological Modelling 141(1-3), 105-123.

AC C

EP

TE D

M AN U

SC

RI PT

626 627 628

29

ACCEPTED MANUSCRIPT Omlin, M. and Reichert, P. (1999) A comparison of techniques for the estimation of model prediction uncertainty. Ecological Modelling 115(1), 45-59.

668 669

Omlin, M., Reichert, P. and Forster, R. (2001b) Biogeochemical model of Lake Zurich: model equations and results. Ecological Modelling 141(1-3), 77-103.

670 671

Pano, A. and Middlebrooks, E.J. (1982) Ammonia Nitrogen Removal in Facultative Wastewater Stabilization Ponds. Journal Water Pollution Control Federation 54(4), 344-351.

672 673 674

Pearson, H.W., Avery, S.T., Mills, S.W., Njaggah, P. and Odiambo, P. (1996) Performance of the phase II Dandora waste stabilisation ponds the largest in Africa: The case for anaerobic ponds. Water Science and Technology 33(7), 91-98.

675 676 677

Privoznik, K.G., Daniel, K.J. and Incropera, F.P. (1978) Absorption, Extinction and Phase Function Measurements for Algal Suspensions of Chlorella-Pyrenoidosa. Journal of Quantitative Spectroscopy & Radiative Transfer 20(4), 345-352.

678 679

Reed, S.C. (1985) Nitrogen Removal in Wastewater Stabilization Ponds. Journal (Water Pollution Control Federation) 57(1), 39-45.

680 681

Reichert, P. (1994) AQUASIM – A tool for simulation and data analysis of aquatic systems. Water Science and Technology 30(2), 21-30.

682 683

Reichert, P. (2005) UNCSIM a computer programme for statistical inference and sensitivity, identifiability, and uncertainty analysis.

684 685

Reichert, P., Borchardt, D., Henze, M., Rauch, W., Shanahan, P., Somlyody, L. and Vanrolleghem, P.A. (2001) River Water Quality Model, IWA Publishing.

686 687

Reichert, P. and Omlin, M. (1997) On the usefulness of overparameterized ecological models. Ecological Modelling 95(2-3), 289-299.

688 689

Reichert, P. and Vanrolleghem, P. (2001) Identifiability and uncertainty analysis of the River Water Quality Model No. 1 (RWQM1). Water Science and Technology 43(7), 329-338.

690 691

Rodi, W. and Research, I.A.f.H. (1984) Turbulence models and their application in hydraulics: a state of the art review, International Association for Hydraulic Research.

692 693

Sah, L., Rousseau, D.P. and Hooijmans, C.M. (2012) Numerical modelling of waste stabilization ponds: where do we stand? Water, Air, & Soil Pollution 223(6), 3155-3171.

694 695

Sah, L., Rousseau, D.P.L., Hooijmans, C.M. and Lens, P.N.L. (2011) 3D model for a secondary facultative pond. Ecological Modelling 222(9), 1592-1603.

696

Saltelli, A. (2002) Sensitivity analysis for importance assessment. Risk Analysis 22(3), 579-590.

697 698

Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M. and Tarantola, S. (2008) Global Sensitivity Analysis: The Primer, Wiley.

699

Shilton, A. (2005) Pond treatment technology, IWA publishing.

700 701

Snee, R.D. (1977) Validation of Regression-Models - Methods and Examples. Technometrics 19(4), 415-428.

702

Spriet, J. (1985) Structure characterization-an overview. IFAC Proceedings Volumes 18(5), 749-756.

AC C

EP

TE D

M AN U

SC

RI PT

666 667

30

ACCEPTED MANUSCRIPT Ulrich, M.M., Imboden, D.M. and Schwarzenbach, R.P. (1995) Masas - a User-Friendly Simulation Tool for Modeling the Fate of Anthropogenic Substances in Lakes. Environmental Software 10(3), 177-198.

706 707

Valero, M.A.C. and Mara, D.D. (2007) Nitrogen removal via ammonia volatilization in maturation ponds. Water Science and Technology 55(11), 87-92.

708 709 710

Van Veldhuizen, H.M., Van Loosdrecht, M.C.M. and Heijnen, J.J. (1999) Modelling biological phosphorus and nitrogen removal in a full scale activated sludge process. Water Research 33(16), 3459-3468.

711 712

Vanrolleghem, P.A., Vandaele, M. and Dochain, D. (1995) Practical Identifiability of a Biokinetic Model of Activated-Sludge Respiration. Water Research 29(11), 2561-2570.

713 714 715

Weijers, S.R. and Vanrolleghem, P.A. (1997) A procedure for selecting best identifiable parameters in calibrating activated sludge model no.1 to full-scale plant data. Water Science and Technology 36(5), 69-79.

716 717 718

Yun, Y.S. and Park, J.M. (2001) Attenuation of monochromatic and polychromatic lights in Chlorella vulgaris suspensions. Applied Microbiology and Biotechnology 55(6), 765-770.

AC C

EP

TE D

M AN U

SC

RI PT

703 704 705

31

ACCEPTED MANUSCRIPT

Table 1 Overview of microbial competition for substrates (S) and dependencies through formed products (P) in the FP. Soluble COD

S

P S

Acetate

NH4+

P S

NO3-

PO43-

O2

N2

CH4

RI PT

Anaerobic populations Hydrolysis Fermenting bacteria (FB) Acetotrophic methanogens (AM) Hydrogenotrophic methanogens (HM) Anoxic and aerobic populations Anoxic heterotrophs (HB) Aerobic heterotrophs (HB) Autotrophs (AB) Algae growth on NH4+ Algae growth on NO3Phosphate accumulating organisms (PAO)

Particulate COD

P P

S S

S S

S P

SC

S S

S

EP

TE D

M AN U

S

AC C

H2

P

S S P P S

P P

S

ACCEPTED MANUSCRIPT ) of the top 30 parameters

Rank 1 2 3 4 5 6 7 8 9 10

Rank 11 12 13 14 15 16 17 18 19 20

Parameters Kz pH Tw YPAO θTw µA YA µH µ ALG YH

4.025 2.267 0.444 0.429 0.428 0.292 0.252 0.196 0.193 0.191

Parameters kz YPHA YFB fp1 iPSF bALG YHM Kh

0.114 0.066 0.059 0.038 0.037 0.036 0.032 0.031 0.021 0.018

Rank 21 22 23 24 25 26 27 28 29 30

Parameters YPO4 YAM bH bHM

Ta iPXB nP iNSF

0.018 0.016 0.013 0.012 0.011 0.010 0.006 0.006 0.006 0.005

RI PT

Table 2. Importance rankings (

AC C

EP

TE D

M AN U

SC

Physical-related parameters: Kz: eddy vertical diffusivity; T w: water temperature; θTw: temperature coefficient; Ta: air temperature. Anaerobic-related parameters: YFB: Yield of FB; YHM: Yield of HM; Kh: Hydrolysis rate; YAM: Yield of AM; bHM: decay rate of HM. Heterotrophs-related parameters: µ H: Maximum growth of heterotrophs; YH: Yield of heterotrophs; bH: decay rate of heterotrophs; : Fermentable substrate half saturation of heterotrophs; : Oxygen half saturation of heterotrophs. Autotrophs-related parameters: µ A: Maximum growth of autotrophs; YA: Yield of autotrophs; : Oxygen half saturation of autotrophs. Algae-related parameters: pH; µ ALG: Maximum growth of algae; : Respiration rate of algae; fp1: Fraction of XI formed during the decay of algae; bALG: Decay rate of algae. PAO-related parameters: YPAO: Yield of PAO; YPHA: PHA requirement for PP storage; iPSF: fraction of P in XS; YPO4: PP requirement (SPO4 release) for PHA storage; iPXB: Fraction of P in bacteria; nP: Correction factor for µ PAO under anoxic conditions

ACCEPTED MANUSCRIPT Table 3. Collinearity index ( ) and determinant measure ( ) of selected parameters subsets. Based on their function, the parameters are categorized into five process groups. The subset in bold is selected for calibration based on the values of and . The description of the parameters can be found in the footnote of Table 2.

Aerobic & Anoxic

PAO

AC C

EP

Combination

Parameters Kz, Tw, kz Kz, Tw, θTw Kz, Tw, θTw, kz YFB, YHM YFB, Kh YFB, YHM, Kh pH, µ ALG, fp1 pH, µ ALG, bALG pH, µ ALG, pH, µ ALG, fp1, bALG pH, µ ALG, , fp1 µ A, YA, µ H µ A, YA, YH µ A, YA, µ A, YA, µ H, YH µ A, YA, µ H, YPAO, YPHA YPAO, iPSF YPAO, YPHA, iPSF Kz, Tw, YFB, pH, µALG, µA, Kh, YPAO, YPHA, YHM, YA, YH Kz, Tw, YFB, pH, µ ALG, µ A, Kh, YPAO, YPHA, bALG, YA, YH Kz, Tw, YFB, pH, µ ALG, µ A, Kh, YPAO, YPHA, µ H, YH, iPSF Kz, Tw, YFB, pH, µ ALG, µ A, Kh, YPAO, YPHA, fp1, µ H, YA Kz, Tw, YFB, pH, µ ALG, µ A, Kh, YPAO, YPHA, bALG, µ H, kz Kz, Tw, YFB, pH, µ ALG, µ A, Kh, YPAO, YPHA, fp1, µ H, YH, kz

RI PT

Algal activity

14.74 11.21 7.91 0.72 0.67 0.62 6.31 5.55 5.53 3.49 3.15 3.69 3.66 3.29 3.86 1.80 2.42 3.02 1.69 3.47 3.37 3.45 3.19 2.51 2.75

SC

Anaerobic

1.24 12.12 12.71 1.13 1.01 1.13 1.20 1.72 3.94 1.84 6.32 6.78 6.80 7.32 6.80 7.34 3.96 1.01 3.97 8.16 8.31 13.43 18.73 29.04 29.06

M AN U

Physical

Set size 3 3 4 2 2 3 3 3 3 4 4 3 3 3 4 4 2 2 3 12 12 12 12 12 13

TE D

Group

ACCEPTED MANUSCRIPT Table 4. Original and calibrated values of the selected parameters.

Original value 0.05 3 2 2 7.8 18.6 0.24 0.053 0.63 0.02 0.63 0.2

Calibrated value 0.07 1.53 1.64 2.98 9.09 25.1 0.12 0.05 0.60 0.02 0.66 0.21

M AN U TE D EP AC C

Change (%) 49.3 -49.0 -17.8 49.1 16.6 34.8 -49.1 -4.6 -4.8 -4.9 4.7 4.9

RI PT

Unit m2.d-1 d-1 d-1 d-1 o C g COD.g-1 N g COD.g-1 COD g COD.g-1 COD g COD.g-1 COD g COD.g-1 COD g COD.g-1 COD

SC

Parameter Kz Kh µA µ ALG pH Tw YA YFB YH YHM YPAO YPHA

ACCEPTED MANUSCRIPT Table 5. Weighted sum of squares (WSS) and mean absolute percentage error (MAPE) for initial and end values (after calibration). MAPE formula can be found in Supplementary Material S8. WSSini 155.9 31.2 6.9 7.0

WSSend 90.1 33.2 3.1 7.4

Change (%) 42.2 -6.3 54.7 -5.6

MAPE (%) 65.9 27.3 14.9 10.7

MAPE (%) 48.9 33.4 9.0 10.2

Change (%) 25.8 -22.2 39.6 4.6

AC C

EP

TE D

M AN U

SC

RI PT

Parameters BOD COD TN TP

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

ACCEPTED MANUSCRIPT

An all-encompassing mechanistic model of a facultative pond is developed

•

A sensitivity analysis reveals the most influential parameters on pond behavior

•

Poor identifiability is caused mainly by the depiction of algal light dependence

•

Model calibration exposes the distinctive characteristics of high-altitude ponds

•

Avoidable uncertainty can be reduced via a systematic parameter estimation

AC C

EP

TE D

M AN U

SC

RI PT

•