Development of a mathematical model for the anaerobic digestion of antibiotic-contaminated wastewater

Accepted Manuscript Title: Development of a Mathematical Model for the Anaerobic Digestion of Antibiotic-Contaminated Wastewater Authors: Rafael Frederico Fonseca, Guilherme Henrique Duarte de Oliveira, Marcelo Zaiat PII: DOI: Reference:

S0263-8762(18)30188-6 https://doi.org/10.1016/j.cherd.2018.04.014 CHERD 3131

To appear in: Received date: Revised date: Accepted date:

24-11-2017 4-4-2018 10-4-2018

Please cite this article as: Fonseca, Rafael Frederico, de Oliveira, Guilherme Henrique Duarte, Zaiat, Marcelo, Development of a Mathematical Model for the Anaerobic Digestion of Antibiotic-Contaminated Wastewater.Chemical Engineering Research and Design https://doi.org/10.1016/j.cherd.2018.04.014 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.

Development of a Mathematical Model for the Anaerobic Digestion of AntibioticContaminated Wastewater

IP T

Rafael Frederico Fonseca*; Guilherme Henrique Duarte de Oliveira; Marcelo Zaiat

Biological Processes Laboratory, Center for Research, Development and Innovation in

SC R

Environmental Engineering, São Carlos School of Engineering (EESC), University of

São Paulo (USP), Engenharia Ambiental - Bloco 4-F, Av. João Dagnone, 1100 - Santa

U

Angelina, 13.563-120, São Carlos, SP, Brazil

M

A

N

*Corresponding author: Phone: 55 16 98136 7853 [email protected]

A

CC E

PT

ED

Graphical abstract

Highlights A mathematical model was developed to describe sulfamethazine degradation



The model was divided into two major stages: acids formation and consumption



Three hypotheses for sulfamethazine degradation were considered



A long-term quantification of the influent variations effects was developed



8 µg of SMZ had a similar impact on the process as 2000 mg of filtered COD

SC R

IP T



Abstract. Anaerobic digestion has been investigated as a potential method for treating antibiotic-contaminated livestock wastewaters. Antibiotic removal is mainly associated

U

with biodegradation and sludge adsorption. In environmental concentrations, i.e., from

N

ng.L-1 to a few hundred µg.L-1, cometabolism is the most likely biodegradation pathway.

A

The overall performance of anaerobic processes may be affected by the hydraulic

M

retention time, and these processes are strongly related to the physical characteristics of

ED

the reactor and variations in influent chemical composition. The effects of these factors can be better understood using a mathematical model. Therefore, this paper aimed to

PT

develop a model to describe an anaerobic process to treat sulfamethazine (SMZ), which was divided into two stages of microorganism growth and substrate consumption. In

CC E

addition, three hypotheses regarding sulfamethazine degradation, including substrate cometabolism related to both stages and an apparent enzymatic reaction, were evaluated.

A

A long-term kinetics structure was added to the model to simulate the process of adaptation to each new operational condition. The results showed that sudden increases in chemical oxygen demand (COD) and hydraulic retention time (HRT) had the most significant negative impact on process performance. In addition, a sudden variation of 8 µg of SMZ had a similar impact on the process as did 1000 mg of filtered COD. Of the degradation hypotheses, the hypothesis related to organic acid consumption was more

accurate than that related to hydrolysis; however, neither could account for the response to variations in HRT. The enzymatic approach resulted in a considerably more accurate representation of the influent flow rate variations than did the cometabolic hypotheses.

IP T

Keywords: Anaerobic digestion, Sulfamethazine degradation, Mathematical modeling,

SC R

Long-term modeling

1. Introduction

U

Antibiotics are among the most successful classes of drugs used to improve health

N

conditions for both humans and animals. Among these drugs, sulfonamides are the most

A

exploited synthetic antibiotic class, and they have been extensively applied to enhance

M

animal feeding and used as growth promoters (Hruska & Franek, 2012). Once consumed, antibiotics are generally transformed by body tissues, where they may be metabolized and

ED

converted to other products that are excreted via urine and feces. However, depending on the substance, approximately 15% to 60% of the original form may be excreted (Perez et

PT

al., 2005). The spread of antimicrobial substances in the environment is a cause for

CC E

increasing concern because of ecotoxicological effects and the possible development of antimicrobial resistance (Baran et al., 2011). Recalcitrant synthetic antibiotics are of interest because their resistance to degradation leads to prolonged exposure times and

A

increased dispersal from their original emission sources. With the development of analytical methods for detecting antibiotics, several

classes have been detected in soil particles, landfills, groundwater, and municipal wastewater treatment effluents (Shelver et al., 2010; Wang et al., 2014; Yin et al., 2014). The presence and concentration of these compounds in the environment also depends on

the degradation capability of wastewater treatment plants (Michael et al., 2013). Most recent studies that have examined pharmaceutical degradation during wastewater treatment were performed on activated sludge systems, due to their widespread application in domestic sewage treatment (Perez et al., 2005). Biotransformation of sulfonamides in aerobic wastewater treatment systems are highly variable, ranging from

IP T

negligible to nearly complete removal (Larcher & Yargeau, 2012; Onesios et al., 2009). These variable results arise from the fact that wastewater treatment systems are primarily

SC R

designed to remove organic carbon and nutrient loads. Therefore, the removal of

pharmaceuticals in these systems occurs as a side effect of this biological activity, and the operational conditions selected to treat the bulk of the organic carbon do not

N

U

necessarily result in reproducible antimicrobial removal efficiencies.

A

Fewer studies have focused on sulfonamide biodegradation in anaerobic treatment

M

systems relative to aerobic systems. Some studies have shown that sulfonamide sulfamethoxazole was completely removed, but sulfamethazine degradation was not

ED

observed (Carballa et al., 2007; Mohring et al., 2009), which indicates that the degradation mechanisms involved are crucial to the performance of the process.

PT

Moreover, antibiotic removal is mainly associated with biodegradation and sludge

CC E

adsorption. Furthermore, in environmental concentrations, i.e., from ng.L-1 to a few hundred µg.L-1, cometabolism is the most likely biodegradation pathway for environmental micropollutants such as the veterinary antibiotic sulfamethazine, based on

A

energy support substrates (Oliveira et al., 2016; Pomies et al., 2013). The physical characteristics of an anaerobic reactor may impact the overall

performance of the degradation process because mass transfer depends on the concentration and characteristics of the biomass (Saravanan & Sreekrishnan, 2006). Moreover, the removal of micropollutants from wastewaters depends on the

physicochemical properties of the substance, and the removal efficiency varies with the reactor’s operational conditions (Pomies et al., 2013). In addition, sudden variations in hydraulic retention time (HRT) or organic loading rates (OLRs) may shift the metabolic pathways between volatile fatty acid producers and consumers, thereby inhibiting acid degradation and influencing the performance of the process (Leitao et al., 2006).

IP T

Mathematical models can be used to understand and quantify how these factors

affect the process. Despite the large number of models in the literature, few models are

SC R

suitable for processes that address antibiotic treatment via anaerobic digestion (AD).

Thus, the aim of this report was to develop a model structure to describe chemical oxygen

U

demand (COD) consumption and volatile fatty acid (VFA) formation and consumption

N

during an anaerobic process treating swine wastewater containing sulfamethazine (SMZ)

A

and to evaluate removal mechanisms under several distinct OLRs and HRTs in a

ED

M

horizontal anaerobic immobilized biomass (HAIB) reactor.

PT

2. Experimental Procedures

CC E

The experiments were conducted in the HAIB reactor previously described by Oliveira et al. (2017). The reactor consists of an acrylic tube with an inner diameter of 5 cm, a length of 100 cm, and a total volume of 1810 mL. The reactor was filled with 0.5-

A

cm cubes of polyurethane foam, which resulted in a bed porosity of 0.62 and a useful volume of 1022 mL. Four equally spaced sampling devices were used to measure the COD, VFA, SMZ and biomass along the length of the reactor. Synthetic wastewater was used to evaluate the process efficiency under ten distinct conditions and during the acclimatization period (Oliveira, 2016). The acclimatization period was divided into two

phases: the first was acclimatization to the OLR range of this study, and the second was acclimatization to the presence of SMZ in the influent (experiments E1 and E2 in Table 1, respectively). The average influent CODt (particulate and soluble), the filtered CODf (soluble) and the SMZ concentrations of each of the operational phases are also shown in Table 1, along with their respective HRT and duration. Notably, in E10, the VFA

IP T

measurement for the first sampling device is missing; however, the other VFA measurements were collected twice.

SC R

3. Mathematical Modeling

The anaerobic process was divided into two stages: 1) hydrolysis/acid formation

U

and 2) acid consumption. SMZ degradation was evaluated during both stages. The model

N

was implemented in Matlab™ and simulated using the stiff solver ODE15s. Among the

A

available solvers, this solver chosen due to its faster response time for each simulation.

M

The hardware used was an Intel I5 4590 running at 3.5 GHz.

ED

3.1. Growth and Substrate Consumption Modeling Several mechanisms were evaluated to match the biomass concentration, total

PT

COD consumption, and acid formation/consumption. Certain assumptions were made to

CC E

achieve the modeling objective, as described below. 1. The process can be represented in two stages: hydrolysis/acidogenesis and acid consumption, and one microbial population is associated with each stage.

A

2. The biomass concentration varies along the reactor. 3. The biomass grows attached to the inert support in overlapping layers, which implies that at the lower layers, the microbes are more strongly attached to the support and consequently will present a larger sludge retention time.

4. Because of bed porosity, the reactional volume varies as a function of the biomass concentration. The HAIB reactor modeling approach was a series of N continuously stirred tank reactors; it was considered to be a plug flow, which was split into 50 reactors (de Nardi et al., 1999). Because the foam matrix porosity was 0.62 and the density of the

IP T

polyurethane foam was 28 kg. m−3 , the reaction volume (Vr) without any biomass was

−1 considered to be 13.57 mL. g Foam . The sludge density was considered to be 1.0 g. mL−1,

SC R

and the reactional volume fraction (η) during the experiment varied along the reactor as a function of the biomass according to Equation (1) as follows: 𝑋

(1)

𝑋 𝑉𝑟

U

𝜂 = 1−𝜌

N

The first stage of the model corresponds to the hydrolysis of the complex substrate and

A

the consequent formation of volatile acids, which are consumed during the second stage.

M

For both stages, the measured values are in units of COD. The influent concentration was

ED

measured as the total COD (CODt) (particulate and filtered) and filtered COD only (CODf). The microbial growth was divided into two populations: a hydrolytic/acidogenic

PT

and an acid-consuming group. The main restriction on growth in both groups was spatial; therefore, the Contois function (Contois, 1959; Wang & Witarsa, 2016) presented the best

CC E

adjustment for substrate degradation and the measured biomass. Equations (2) and (3) represent the CODf consumption and acid formation/consumption stages, respectively,

A

and Equations (4) and (5) determine the substrate consumption rate for each stage. Equation (6) represents the liquid-phase mass transfer coefficient of the HRT, in which vs is the liquid superficial velocity, as defined in Equation (8). 𝑑𝐶𝑂𝐷𝑘 𝑑𝑡

= (𝐶𝑂𝐷𝑘−1 − 𝐶𝑂𝐷𝑘 )

𝐷𝐻 𝜂

− 𝐹1 𝑘1 𝑆1𝑘 𝑋ℎ𝑘

(2)

𝑑𝑉𝐴 𝑘 𝑑𝑡

𝑘1 = 𝑘2 =

= (𝑉𝐴 𝑘−1 − 𝑉𝐴 𝑘 )

𝐷𝐻 𝜂

𝐾1𝑚 𝐾𝑆ℎ (𝑋ℎ𝑘 +𝑋𝑎𝑘 )+𝑆1𝑘 𝐾2𝑚 𝐾𝑆𝑎 (𝑋ℎ𝑘 +𝑋𝑎𝑘 )+𝑉𝐴𝑘

+ 𝐹1 𝑘1 𝑆1𝑘 𝑋ℎ𝑘 − 𝐹2 𝑘2 𝑉𝐴𝑘 𝑋𝑎𝑘

(3)

𝑘𝑠

(4)

𝑘𝑠

(5)

𝑘𝑠 = 0.033𝑒 0.0217𝑣𝑠

IP T

(6)

where S1k = CODk − VAk − CODe , with CODk (mg CODf . L−1 ) representing the measured

SC R

CODf and CODe representing the inert part of the influent CODf, which was estimated to

be approximately 7%. Therefore, S1k (mg CODf . L−1 ) was the hydrolyzable mass of the

U

−1 CODf. The function parameter K1,2m (mg CODf . g F . L−1 g −1 SSV cm ) and the parameter

N

K Sh,a (mg CODf . g F . g −1 SSV ) represent the maximum rates for hydrolysis and acid

A

consumption/formation and the half-saturation biomass for the Contois kinetics,

M

respectively. The kinetic parameter K1,2m is the substrate consumption constant, which was adjusted for each experiment, that in association with k s (cm. h−1 ), which is the

ED

liquid-phase mass transfer expression; determine the substrate consumption rate. According to (Sarti et al., 2001), ks affects the overall reaction rate effectiveness and thus

PT

was incorporated in the modeling to represent the HRT effects. The function was

CC E

calculated in a similar HAIB but with a bed porosity of 0.40. Note that K1,2m is also included in the OLR function, which is defined in the Long-Term Quantification of OLR and HRT Effects section, and F1,2 is the non-dimensional function that represents biomass

A

acclimatization/mutation/specialization (Ghasimi et al., 2015) in the process OLR, which ranges from zero to one, as shown in Equations (18) and (19). DH (h−1 ) is the dilution rate for each tank and is calculated as the ratio of the liquid superficial velocity vs (cm. h−1 ), which is defined in Equation (7). The DH is described by Equation (8), where

L=100 cm is the reactor length, and N=50 is the number of tanks. The HRT varied, as shown in Table 1. 𝐿

𝑣𝑠 = 𝐻𝑅𝑇

(7)

𝑣

𝑠 𝐷𝐻 = 𝐿/𝑁

(8)

IP T

The biomass concentration for each tank was defined as follows in Equations (9) and (10):

𝑑𝑡

𝑑𝑋𝑎 𝑘

𝐷𝑠

= Ya 𝐹2 𝑘2 𝑉𝐴𝑘 𝑋𝑎𝑘 − (𝐾𝑑 +

𝐷𝑠

𝜂

𝜂

𝑄 2

(𝑄 ) ) (𝑋ℎ𝑘 − 𝑋ℎ 𝑟𝑒𝑠 ) 0

𝑄 2

( ) ) (𝑋𝑎𝑘 − 𝑋𝑎 𝑟𝑒𝑠 ) 𝑄0

(9)

(10)

U

𝑑𝑡

= Yh 𝐹1 𝑘1 𝑆1𝑘 𝑋ℎ𝑘 − (𝐾𝑑 +

SC R

𝑑𝑋ℎ 𝑘

N

−1 where Yh,a (g SSV . L. g F−1 . mg −1 COD ) is the growth yield coefficient, and Ds (h ) is the

A

2 dilution rate for all the experiments and is associated with the term η−1 (Q. Q−1 0 ) (Horn

M

et al., 2003), which represents the effects of the shear stress force on the sludge retention time caused by the HRT and the flow variations in each tank caused by the biomass

ED

concentration. In Equations (9) and (10), the coefficient Q represents the influent volumetric flow, e.g., Q = V⁄HRT (L. h−1 ), and Q0 is the volumetric flow for an HRT of

PT

24 h. Xh res and Xa res are the biomass amounts expected to be inactive during the removal

CC E

process in the near absence of energy substrate, and they were assumed to be X0 Yh ⁄(Yh + Ya ) and X 0 Ya ⁄(Yh + Ya ) for stages 1 and 2, respectively, where X0 is the

A

measured initial biomass content. Lastly, K d (h−1 ) is the biomass death rate.

3.2. Sulfamethazine Degradation Three hypotheses were evaluated for SMZ degradation. The first two hypotheses were related to cometabolism, as proposed by (Criddle, 1993), and these hypotheses are

associated with both stages. Hydrolysis-associated cometabolism (first hypothesis) and acid-consumption-associated cometabolism (second hypothesis) are shown in Equations (11) and (12), respectively, as follows:

𝑑𝑡 𝑑𝑆𝑀𝑍𝑘 𝑑𝑡

= (𝑆𝑀𝑍𝑘−1 − 𝑆𝑀𝑍𝑘 )

𝐷𝐻

= (𝑆𝑀𝑍𝑘−1 − 𝑆𝑀𝑍𝑘 )

𝐷𝐻

𝜂

𝜂

− 𝑌𝑧1 𝑘1 𝑆1𝑘 𝑋ℎ𝑘 𝑆𝑀𝑍𝑘

(11)

− 𝑌𝑧2 𝑘2 𝑉𝐴𝑘 𝑋𝑎𝑘 𝑆𝑀𝑍𝑘

(12)

IP T

𝑑𝑆𝑀𝑍𝑘

SC R

where YCh,a (𝑚𝑔𝑆𝑚𝑧 ⁄𝑚𝑔𝐶𝑂𝐷 ) is the SMZ cometabolic transformation capacity.

The third hypothesis is shown in Equations (13) and (14) and is expressed as the apparent enzymatic activity degradation. This hypothesis proposes that enzymes are produced

U

during the hydrolysis stage and secreted to the bulk liquid, where the reactions occur. The

N

acid-consuming stage was not considered for this hypothesis because the associated

follows:

𝑑𝑡 𝑑𝐸𝑀𝑍𝑘

𝐷𝐻

𝜂

− 𝑉𝐸 𝐸𝑀𝑍𝑘 𝑆𝑀𝑍𝑘

(13)

= (𝐸𝑀𝑍𝑘−1 − 𝐸𝑀𝑍𝑘 )

𝐷𝐻

𝜂

− 𝑌𝐸 𝑘1 𝑆1𝑘 𝑋ℎ𝑘 − 𝐾𝐷𝑒 𝐸𝑀𝑍𝑘

(14)

PT

𝑑𝑡

= (𝑆𝑀𝑍𝑘−1 − 𝑆𝑀𝑍𝑘 )

ED

𝑑𝑆𝑀𝑍𝑘

M

A

enzymes are expected to be intracellular (Batstone et al., 2002). The equations are as

CC E

where EMZk (U. L−1 ) is the apparent enzymatic activity, VE (L. U −1 . h−1 ) is the maximum SMZ degradation rate constant, YE (U. mg −1 COD ) is the apparent enzymatic activity

A

formation yield constant, and K De (h−1 ) is the enzymatic activity degradation rate.

3.3. Long-Term Quantification of OLR and HRT Effects The effects of operational condition variations modeling were considered in three steps. The first two focused on parameter adjustments for the proposed model; in the third step,

sensitivity and collinearity analyses were performed to avoid overparameterization and exclude unnecessary parameters. The first step focused on the amplitude of the parameters due to the influent OLR. Equations (15) through (17) express these variations for k1 , k 2 , Y𝑧 1,2 and YE as calculated. K1M and K 2M are the maximum consumption rates for their

respective functions, and YSM represents the maximum transformation capacity

IP T

coefficient for Yz1,2 and YE , which is represented by Yz in Equation (20). K 𝑠1 , K 𝑠2 and K zi

are the half-saturation constants, and the index i in K zi refers to the Y indices. For the

SC R

cometabolic hypothesis, the parameter units are (μg Smz . mg −1 COD ), and for the apparent f

enzymatic activity, the parameter units are (U. mg −1 COD ).

𝑠1 +𝐶𝑂𝐷0

𝑌𝑠 = 𝐾

𝐶𝑂𝐷0 𝑠2 +𝐶𝑂𝐷0

)

M

𝐾2𝑚 = 𝐾2𝑀 (𝐾

)

(15)

N

𝐶𝑂𝐷0

A

𝐾1𝑚 = 𝐾1𝑀 (𝐾

U

f

𝑌𝑆𝑀

(17)

ED

𝑧𝑖 +𝐶𝑂𝐷0

(16)

The second step was performed to adapt the populations to a new operational condition,

PT

which is expressed by a logistic equation and subordinated to changes in the OLR and

CC E

HRT by three sensitivity functions, as shown in Equations (18) to (20). In this sense, the growth parameters in Equations (4) and (5) are now a function of Equations (15) to (17) with the population adaptation dynamics in Equations (18) to (20). Here, F1 and F2

A

represent the overall activity of hydrolytic/acid formation bacteria and acid-consuming bacteria (Dalsenter et al., 2005), and Yzi represents the overall transformation capacity of the SMZ-associated bacteria, with the index i representing each of the degradation hypotheses. The values are expressed as a percentage. The parameter fvz refers to the decrease rate of Yzi , which is lower than the increase rate. The equations are as follows:

𝑑𝐹1 𝑑𝑡 𝑑𝐹2 𝑑𝑡

𝐹

1 = 𝜈1 𝐹1 (1 − 𝐹1 ) − (Δ𝐶𝑂𝐷0 𝑘1 + Δ𝑆𝑀𝑍0 𝑘1 + Δ𝑣𝑠𝑘1 ) 100

𝐹

2 = 𝜈2 𝐹2 (1 − 𝐹2 ) − (Δ𝐶𝑂𝐷0 𝑘2 + Δ𝑆𝑀𝑍0 𝑘2 + Δ𝑣𝑠𝑘2 ) 100

𝑌

𝑑𝑌𝑧 𝑑𝑡

(18)

(19) 𝑌

𝑑𝑌

𝑧 𝑓𝑣𝑧 𝜈𝑧 𝑌𝑧 (1 − 𝑌𝑧 ) − (Δ𝐶𝑂𝐷0 𝑌𝑠 + Δ𝑆𝑀𝑍0 𝑌𝑠 + Δ𝑣𝑠𝑌𝑠 ) 100 , 𝑑𝑡𝑧 < 0 𝑠 ={ 𝑌𝑧 𝑌𝑧 𝑑𝑌𝑧 𝜈𝑧 𝑌𝑧 (1 − 𝑌 ) − (Δ𝐶𝑂𝐷0 𝑌𝑠 + Δ𝑆𝑀𝑍0 𝑌𝑠 + Δ𝑣𝑠𝑌𝑠 ) 100 , 𝑑𝑡 ≥ 0

(20)

IP T

𝑠

where the sensitivity functions are expressed in Equations (21) to (23). ΔCOD0 i represents

SC R

the OLR shock sensitivity to the hydrolysis, acid-consuming and SMZ degradation

stages, which are indicated by the index i. Note that in all three cases, CODt was considered for the OLR shock loads, as shown in Table 1. ΔSMZ0 i and Δvsi reflect the

U

process sensitivity for SMZ and HRT shocks, respectively. In Equations (21) and (22),

N

the f1,2,z value represents the sensitivity fraction when a negative variation of any effluent

A

concentration occurs, i.e., the present value is lower than prior values. The K e1,2,z , p1,2 and

M

pSmz values determine the sensitivity of the populations to OLR changes, and higher

ED

values correspond to lower sensitivity to small changes in OLR. K Sz i and fzi determine the sensitivity to SMZ variations and K vi and fvi to HRT. Δ𝐶𝑂𝐷0 (𝑡)𝑝 𝐾𝑒𝑖

Δ𝐶𝑂𝐷0 (𝑡)𝑝

CC E

𝐾𝑒𝑖

A

Δ𝑆𝑀𝑍0 𝑖 = {

Δ𝑣𝑠𝑖 = {

𝑓𝑧

𝑖

, Δ𝐶𝑂𝐷0 (𝑡) ≥ 0

Δ𝑆𝑀𝑍0 (𝑡)𝑝𝑠𝑚𝑧

𝑖

𝐾𝑆 𝑧

𝑖 Δ𝑆𝑀𝑍0 (𝑡)𝑝𝑠𝑚𝑧

𝐾𝑆 𝑧

𝑓𝑣

, Δ𝐶𝑂𝐷0 (𝑡) < 0

PT

Δ𝐶𝑂𝐷0 𝑖 = {

𝑓𝑐 𝑖

Δ𝑣𝑠 𝐾𝑣𝑖

Δ𝑣𝑠 𝐾𝑣𝑖

𝑖

(21)

, Δ𝑆𝑀𝑍0 < 0

, Δ𝑆𝑀𝑍0 ≥ 0

(22)

, Δ𝑣𝑠 (𝑡) < 0

, Δ𝑣𝑠 (𝑡) ≥ 0

(23)

Because of certain observations of the effluent response to OLR shocks in the influent, Equations (24) to (26) were added to avoid sudden responses in the process simulations. Thus, the effects of HRT and OLR variations in the overall performance of the process were quantified by the time (in days) required for the process to be affected by those changes. The parameters n and m represent the number of days for the effluent to react to

IP T

the influent variations. Similar behavior was studied by (Amorim et al., 2005), in which the time required for the process to return to the previous equilibrium was evaluated when

SC R

applying different OLR shock load impulses. 1

Δ𝐶𝑂𝐷0 (𝑡) = 𝑛 ∑𝑛𝑖=0 (𝐶𝑂𝐷0 (𝑡 − 𝑖) − 𝐶𝑂𝐷0 (𝑡 − 𝑖 − 1))

(24)

(25) (26)

M

A

N

1

Δ𝑣𝑠 (𝑡) = 𝑚 ∑𝑚 𝑖=0(𝑣𝑠 (𝑡 − 𝑖) − 𝑣𝑠 (𝑡 − 𝑖 − 1))

U

1

Δ𝑆𝑀𝑍0 (𝑡) = 𝑛 ∑𝑛𝑖=0 (𝑆𝑀𝑍0 (𝑡 − 𝑖) − 𝑆𝑀𝑍0 (𝑡 − 𝑖 − 1))

ED

3.4. Parameter Adjustment Procedure

3.4.1. Growth and Substrate Consumption and Sulfamethazine Degradation

PT

Because the work presented here represents a long-term experimental procedure

CC E

and the operational conditions were changed after a certain period, the adjustment procedure was conducted iteratively. The fitting procedure was conducted in three major stages. During the first stage, the structure was defined, and its parameters were adjusted

A

to experiment 7 data because of the amount of data related to the OLR, VFA, and biomass (primarily) along the reactor. During previous experiments, VFAs were not detected during the measurements. The initial conditions of the simulation were the linearly interpolated experimental data of the previous experiment. Once the model structure had

been defined, the next step was to state the initial values for the parameter adjustment procedure. The second stage focused on manually adjusting the parameters for all the experiments, beginning with biomass acclimatization (experiment 1), in which the biomass initial conditions were defined as 𝑋ℎ = 𝑋𝑎 = 0.25 𝑔𝑥 𝑔𝐹−1, and all the other state

IP T

variables were set to zero. The parameters related to biomass growth, i.e., Yh , Ya , K d , Ds , K Sh and K Sa , were considered fixed for all OLR changes. The function

SC R

parameters 𝑘1 and 𝑘2 were sequentially adjusted for each experiment from 1 to 7, and the value of the biomass growth parameters was recalculated at each iteration based on the

U

fitting procedure described in this section. The relationship between the sludge retention

N

time (SRT) and the HRT was defined during this stage.

A

In the third stage, a response surface methodology determined by a central

M

composite design (CCD) technique was associated with an evolutive operation performed to optimize the parameter fitting and suppress any parameters that were not statistically

ED

relevant to the process (Kumar et al., 2011; Singh et al., 2005). This last stage also

PT

followed the iterative idea of the fitting procedure during parameter adjustment. The main difference was that the parameters in this study were adjusted in pairs, e.g., Yh and Ya ; K d

CC E

and Ds ; and K Sh and K Sa . In general, the CCD was conducted with five levels for each pair, although in certain cases, seven levels were required to obtain a more precise assessment of the factors. For the parameters 𝐾𝑑 and 𝐷𝑠 , the main restriction was the

A

variation in biomass between experiments 7 and 11. The cost function was used to maximize the R2 value of the estimations. Once the growth and substrate consumption parameters were fit, the SMZ degradation hypotheses were evaluated. For the Criddle cometabolism (Criddle, 1993),

which is shown in Equations (10) and (11), the YSh and YSa values were simply adjusted to the best R2 for all the experiments. For the enzymatic degradation, the response surface methodology was used to maximize the R2 value. Fitting Procedure

IP T

The fitting procedure was conducted according to the flowchart presented in Figure 1.

SC R

3.4.2. OLR and HRT Changes

First, the parameters of Equations (18) to (20) were simultaneously adjusted with those of Equations (15) and (16) to the calculated values of K1m and K 2m . The values of

N

U

K S1,2 were obtained as the half-saturation constants for carbohydrate hydrolysis presented

A

in the Anaerobic Digestion Model number 1 (ADM1) (Batstone et al., 2002). A manual

M

procedure similar to the fitting procedure previously described, was conducted to adjust the values of ν1,2 , K S1,2 and p1,2 in Equations (21) to (23), and the values of n and m in

ED

Equations (24) to (26). After a satisfactory correlation between K1m and K 2m was found for experiments 2 to 8, the HRT shock parameters K v1 , K v21 , K v22 , and K vs were evaluated

PT

for experiments 9 to 11. An iterative adjustment was performed, and the parameter set of

CC E

Equations (21) to (23) was adjusted based on the calculated K1m and K 2m values; in addition, n and m were adjusted based on the effluent data. The K1M and K 2M values were estimated based on the calculated k1 and k 2 values. After each iteration, a long-term

A

simulation was conducted to correlate the simulated CODf and VFA along the reactor with the measured data of each experiment. This procedure was repeated until a better R2 value could not be found for each experiment. After adjusting the substrate consumption parameters, the SMZ degradation hypotheses were considered. The fitting procedure was analogous to the previously

described procedure except for K sz , which was adjusted based on the previously calculated YSh and YE values. Concerning the second SMZ degradation hypothesis, VE (L. U −1 h−1 ) in Equation (13) was assumed to be 1.0 L. U −1 h−1 because this approach described the apparent enzymatic activity of SMZ degradation, and information

IP T

concerning this activity is not available in the literature.

SC R

3.4.3. Process Identifiability

Once a model has been developed that sufficiently describes the measured data, it is essential to assess whether the model parameters are identifiable by the amount of

U

experimental data. Another critical question is how model predictions are affected by the

N

lack of data (Raue et al., 2009). Therefore, the model identifiability was assessed in two

A

parts. The first concerned growth and substrate consumption and sulfamethazine

M

degradation modeling; the second was related to the long-term quantification of OLR and

ED

HRT effects. In the first step, the proposed model goodness-of-fit was analyzed using the chi-square (𝜒 2 ) method. Thus, considering the number of parameters and the number

PT

experimental measurements, if the calculated chi-square value was lower than a critical value, shown in Table 2, the process was considered identifiable. In the second step, given

CC E

the large number of parameters and measurements, sensitivity and collinearity analyses were performed to evaluate which parameters did not influence to the outputs or if their

A

influence could be compensated for by another parameter in the set (Gabor et al., 2017). The chi-square value was used as measurement of the goodness-of-fit (Vera et al., 1992). The chi-square value was calculated using Equation (27), and the sensitivity and collinearity analyses were conducted based on the guidelines presented by Gabor et al. (2017). Note that for these analyses, the Matlab™ solver ODE113 options set had to be

fixed with a maximum integration step of 0.1 days and both relative and absolute accuracies of 10-7. This solver was used because of the stringent error tolerances required by the problem. Otherwise, the simulation noise would exceed the differences between the parameter variations. In contrast to the procedure proposed by Gabor et al. (2017), the largest identifiable subset was defined manually using the Pearson correlation to identify

IP T

the parameters that should affect the collinearity (Villaverde & Banga, 2014). If a

parameter had a higher correlation index with any other parameter, then it was excluded,

𝜒𝑗2 (𝜃)

=

(𝑦𝑗𝑖 −𝑦̂𝑗𝑖 ) ∑𝑚 𝑖=1 𝜎2

SC R

until the collinearity index was achieved or the removal affected any of the other outputs. 2

(27)

𝑗𝑖𝑆

N

U

where 𝑦𝑗𝑖 are the measured data from operational phase j at the sampling devices i, 𝑦̂𝑗𝑖 are the estimated data from the modeling, 𝜃 is the set of parameters used to adjust the

M

A

model to the experimental observations, and 𝜎𝑗𝑖2 is the sample variance of the substrates COD and SMZ calculated using Equations (28) and (29). 𝜎𝑗𝑖2 was calculated based on the

ED

effluent data, which was normalized and leveled for each sampling point, as shown in Equation (29). The same procedure was conducted for the SMZ data. For VFA, since E2,

PT

E6 and E8 had the same operational conditions, their data were considered in the variance

CC E

calculations. However, because the measured VFA in E6 was zero, its data were discarded; otherwise the VFA variance would have values similar to the measured data. Normalization and leveling procedures were used to estimate the VFA variance during

A

each phase; however, the variance was based on the ratio between the influent CODf of each phase and the 2000 mg O2 . L−1 . d−1 of the respective experimental phases. The equations are as follows: 𝜎𝑗𝑖2𝑒𝑓 = ∑𝑛𝑗=2

(𝑦𝑗𝑖 −𝑦̅𝑗𝑖 ) 𝑛−1

2

(28)

𝜎𝑗𝑖2 =

𝑦𝑗𝑖 𝑦̅𝑗

𝜎𝑗𝑖2𝑒𝑓

(29)

where 𝑦̅𝑗 is the mean of the experimental effluent data for each operational phase, and n=11 is the number of operational phases. The value of j in Equation (28) started at 2

IP T

because of the available experimental data, as shown in Table 2.

In Table 2, Np is the number of parameters, Ns is the number of measured samples, df is

SC R

the degrees of freedom (𝑑𝑓 = 𝑁𝑠 − 𝑁𝑝 ), and CV (95%) represents the chi-square critical

U

value for a 95% confidence interval.

N

4. Results and Discussion

A

In a previous study, Oliveira et al. (2017) evaluated changes in the OLR and HRT

M

and their effects on the process stability and performance for both COD and SMZ

ED

removal, and they concluded that reducing the HRT considerably impacted the SMZ degradation but maintained the COD removal efficiency. Moreover, they found that the

PT

OLR had nearly no effect on the overall process performance. By modeling the data obtained by Oliveira et al. (2017), other hypotheses concerning this process can be

CC E

considered.

4.1. Biomass Growth and Substrate Consumption

A

As discussed by Donoso-Bravo (2011) (Donoso-Bravo et al., 2011), two questions

arise during modeling: first, can the model structure fit the experimental data without overparameterization, and second, does the model structure allow for a unique optimal set of parameters? Concerning the first question, several model structures were tested; however, only the combination of the reactional volume fraction 𝜂 with the Contois

kinetics resulted in an accurate representation of the state variables of the experimental conditions. The Contois function is commonly used to describe hydrolysis/disintegration of particulate wastes (Ramirez et al., 2009), e.g., when biomass attaches to an insoluble particle and grows around it. As a consequence, the outer biomass cellular activity is reduced. However, in the reactor used here, the biomass attaches to the foam matrix and

similar.

With

respect

to

the

second

question,

and

still

IP T

the biological activity of the inner layers of biomass is weakened; thus, the effects are considering

the

SC R

overparameterization problem, during the third stage fitting, K d and Ds and K Sh and K Sa

had to be reevaluated because some of these parameters did not have a significant effect on the modeling accuracy. Therefore, an ambiguity was observed in the optimum

U

adjustment of the biomass when analyzing K d and Ds for experiments 7 and 11 (Table 1).

N

For experiment 7, variations in K d did not have a significant impact; meanwhile, the lower

M

A

the Ds , the higher the R2. In contrast, for experiment 11, the highest R2 values were found with higher Ds values and lower K d values. Another evaluation was conducted but with

ED

a K d near zero. In this situation, only Ds was significant for both experiments. When suppressing the parameter K d , no loss of biomass adjustment was observed for all the

PT

experiments. Thus, the parameter K d was excluded from the modeling.

CC E

Concerning the biomass growth yield coefficients, an optimum value was found for all the experimental conditions. These values are shown in Table 3. For K Sh and K Sa , an optimum was found for experiment 7 at K Sh = 530 and K Sa = 530 mg COD . g F . g −1 SSV .

A

For experiment 11, only the first parameter was significant, and increasing it resulted in better fitting accuracy. Based on these observations, the number of parameters was reduced. K Sh was set equal to K Sa , and a stage three analysis was conducted considering Ds to evaluate the combined effects. An optimum adjustment was found using this analysis, and the results were considered to be the new fixed values for these parameters;

accordingly, the parameters k1 and k 2 were readjusted to fit the experimental observations. During the adjustments for E11, the biomass could only be adjusted if the combined effects of K d and Ds with the liquid mass transfer coefficient were as high as those observed in E9. This observation is discussed in section 4.3.1. The results are shown in Table 3, where the average filtered influent OLR is also described. (COD f ), R2CODt and in

IP T

R2VA are the fitting adjustments of both stages, and χ2CODf and χ2VA are the calculated chi-

SC R

square for each experiment. For experiment 1 of Table 3, no VA data were acquired.

The variations in the parameters k1 and k 2 throughout the experiments were significant, as shown in Table 3. In both cases, these variations were correlated with the

U

process acclimatization and influent concentration. For the hydrolysis stage, the influent

N

concentration was more significant than the acclimatization except in experiment 11, and

A

the opposite was true for the acid-consumption stage. These findings indicate that the first

M

stage is much more robust to process variations than the second because there are more

ED

possible limiting steps in the second stage than in the first (Batstone et al., 2002). Except for E4 and E10, all the chi-square values were within the critical value

PT

limits. For E4, chi-square analysis indicated that hydrolysis stage modeling could not identify the process. This was expected, because as shown in Figure 4c of section 4.3.1,

CC E

the CODf values rose along the reactor, while the first order model effluent decreased, and its value stabilized. Because this was a unique occurrence among the experiments, no

A

further analysis will be conducted. Concerning E10, the lack of the first round of VFA data compromises model identifiability. However, because the other measurements were taken twice, yielding a mean and a variance for each sample, and all the CODf points were collected, the VFA estimation is plausible.

4.2. Sulfamethazine Degradation Hypotheses A comparison between the two direct cometabolism approaches shows that the stage associated with VFA consumption presented a better adjustment than did the hydrolytic stage. This observation is directly linked to the available substrate for both stages because the consumption rate during both stages is related to their substrate

IP T

concentration. In the first stage, the concentration decreased across the reactor length, whereas in the second stage, the concentration was higher near the first sampling port.

SC R

Thus, SMZ degradation was expected to be better represented during the second stage.

However, none of the stages could accurately represent the degradation at a lower HRT.

U

In turn, the enzymatic degradation approach yielded a more accurate reflection of the

N

experiments after the OLR shock of experiment 7, especially for those with lower HRT,

A

as shown in Table 4. The last column of Table 4 is the sum of the chi-square values for

M

all the experiments, and lower values indicate more accurate model predictions. The enzymatic degradation rate was estimated to be 0.40 h−1 . This result balances the faster

ED

response observed in experiments E2 to E6 and the slower responses in experiments E7 to E11. Higher values reduced the contribution of the latter, while lower values reduced

PT

the contribution of the former experiments. Note that the chi-square results suggest that

CC E

acid-consuming cometabolism is the most plausible hypothesis for SMZ degradation, even with the low R2 values in E9 and E10. Notably, the SMZ degradation was considered identifiable in E10 despite the low R2 value. This occurred because the chi-square analysis

A

takes into consideration the process variance, while the R2 calculation only considers the model’s deviations from the experimentally measured means. Thus, the chi-square value is a more robust way to evaluate this process. 4.3. Long-Term Quantification of the OLR, HRT and SMZ Effects 4.3.1. Growth and Substrate Consumption

Table 5 shows the growth and substrate consumption parameters related to the OLR and HRT shock sensitivities. The “Estimated Value” column shows the parameters adjusted to Equations (1) to (26), while the “After S/C Analyses” column shows the parameters adjusted after the removal of non-influential parameters following sensitivity and collinearity analysis, which are shown in Table 6. Note that, for a parameter be

IP T

considered influential to the process, the model must be sensitive to it, and its effects must be unique when compared to those of other parameters (Gabor et al., 2017). In the “After

SC R

S/C Analyses” column of Table 5, the infinity symbol and the three zeros indicates that a

parameter could be removed without any loss to the model goodness-of-fit. Considering the sensitivity and collinearity analyses, the COD and VFA lines in Table 6 are related to

U

the first and the second stages of the model, respectively. In these lines, only the

N

parameters marked with an “S/C” and highlighted in gray were maintained in the model.

A

The collinearity indices for the first and second stages were 15.53 and 14.61, respectively,

M

which were below the critical value of 20.0 proposed by (Gabor et al., 2017).

ED

According to (Ghasimi et al., 2015), reducing or increasing the OLR significantly affects the microbial community genetics. The adaptation/specialization coefficient can

PT

mathematically represent the effects of population changes. Considering the adjusted

CC E

parameters, the adaptation velocity behavior was similar to the genetic specialization reported by Ghasimi et al. (2015). The maximum consumption rate during the first stage was higher than the values shown in Table 3. This difference was related to parameter

A

exclusion. During the collinearity analysis, a linear dependency was observed between the OLR shock load sensitivity parameter and the adaptation rate constant. The former was chosen to be removed because if the latter was removed, no effect would be observed (all the shock recoveries depend on it). In contrast, the second stage community did not

achieve maximum performance, mainly because of the OLR shock loads during the experiments, but also due to SMZ shock sensitivity. It is also important to ensure that the positive and negative HRT shock sensitivity parameters were neither influential to the COD nor to VFA but that they affected the biomass, as shown in Table 6. As shown in Table 3, during E11, the Ds value had to be

IP T

corrected so the process could represent both COD consumption and biomass

concentration along the reactor. However, given a proper value for the positive HRT

SC R

shock sensitivity parameter, the deviation in Ds that should be related to shear stress was compensated for by the biomass recovery. This finding may be related to the liquid mass

U

transfer coefficient from Equation (6), which was reduced from 5.89×10-3 to 1.96×10-3

N

cm.h-1 (Sarti et al., 2001) for an HRT varying from 8 to 24 h. This reduction may have

A

restricted microbial feeding mechanisms, thereby leading to reduced activity. The impacts

K v1,2 in Table 5 and in Figure 3.

M

of positive HRT variations were higher during the first stage, as shown by the parameter

ED

The process response time to OLR, SMZ and HRT variations had a strong

PT

influence on the modeling. These parameters affected the effluent CODf more than the substrate consuming kinetics; these values were measured before the transition to a new

CC E

experiment, as shown in Figure 2. Without the soft transitions, the estimated CODf peaks after E1, E10 and E11 would be twice as high as those observed. If these soft transitions could be considered instantaneous, their behavior would be compensated for by other

A

parameters. However, this compensation would negatively affect the estimates of biomass and substrate consumption kinetics. The Contois inhibition structure was successfully used to represent the spatial restrictions caused by the biomass concentration in the foam matrix. Both half-saturation constants calculated in the ADM1 (Batstone et al., 2002) yielded accurate representations

of the effluent responses to the influent CODf concentration fluctuations, as shown in Figure 2 (a). The simulated and experimental effluent CODs are shown in Figure 2 (b). Note that the differences between the simulated and experimental CODs during E0 (Figure 2 (b)) are due to the assumed inert COD, which is considered to be 7% of the influent CODf (Figure 2 (a)). These differences are also affected by F1 and F2, the initial

IP T

conditions of Equations (18) and (19), which determine the process performance. The

initial conditions were adjusted to 0.30 and 0.015, respectively, for F1 and F2, to match as

SC R

closely as possible the estimated values of k1 and k2 (Table 3, and as shown in Figure 3). Moreover, note that the experimental effluent COD only changes to the assumed inert

U

COD almost 25 days after the OLR reached 2000 (mgO2 . L−1 . d−1 ) during the

N

acclimatization period. The reason for this behavior was not clear. During all the other

A

experiments, the simulated effluent remained within the 95% confidence interval bounds.

M

Note that in E10, the variations in the effluent were of such magnitude that the confidence

ED

intervals could not be shown in Figure 2, which was limited to 500 mgO2 . L−1 . Concerning the effects of OLR changes on the VFA consumption rate, negative

PT

variations had a lesser impact on process performance than did positive variations. This result occurred because although negative variations are related to slow biomass death,

CC E

positive variations are related to microbial adaptation to the new operational conditions, as well as inhibition mechanisms (Hierholtzer & Akunna, 2014). During the model development, the process appeared to be more sensitive to CODt variations than to the

A

filtered COD. The parameters K e2 and p reflect microbial sensitivity to shock. In this sense, lower values for the parameter K e2 correspond to higher shock impacts on the process performance, and higher p values correspond to greater effects associated with

greater OLR changes. In addition, the association between higher K e2 values and higher p values indicates that the microbial community is only sensitive to high OLR changes. The most valuable information about the effect of SMZ variations on the overall process performance was obtained in experiment 2, in which no other significant influent variations were observed. In fact, 8 µg of SMZ shock presented an effect greater than the

IP T

1000 mg CODf shock load for the acid-consuming stage, as observed in the transition from E4 to E5. During the other experiments, the SMZ variations were not high enough to

SC R

disturb the green line in Figure 3. However, the combined effect of small variations affected the second stage overall performance. This result would have been observed in

U

E3 to E7 if the effect were removed; otherwise, the process performance would have

N

become high enough for the model to become unidentifiable based on the chi-square

A

values.

M

Figure 4 compares the experimental data and the simulated data along the reactor length. For the first stage, except for E4, with a goodness-of-fit of 28.75, the other

ED

experiments all showed goodness-of-fit values below 3.97. For the second stage, all

PT

goodness-of-fit values were below 3.43. These results were consistent with those in Table 3 and demonstrate the accuracy of the long-term model in explaining the process behavior

CC E

under several influent compositions. As a consequence, experiments 7 and 11 had R2 values of 0.97 and 0.85, respectively, between the modeled and the measured total biomass (Figure 5). Chi-square values were not calculated because no measure was

A

repeated. In Figure 5, the blue line represents the total estimated biomass in gSSV.gF-1; this value can be compared to the purple ‘x’ marks, which show the measured biomass content for each sampling device. The red and the yellow lines represent the biomass related to the hydrolysis and the acid-consuming stages, respectively. Note that the hydrolysisrelated biomass has a higher concentration at the beginning of the reactor and decreases

along its length, while the acid-consuming related biomass has a low initial concentration, increases near the first sampling device and decreases near the others. These behaviors follow the VFA concentrations shown in Figure 4 and are consistent with the difference between the COD and the VFA (the hydrolyzable part of the substrate). The identifiability

IP T

of this modeling approach is discussed in section 4.4.

SC R

4.3.2. Sulfamethazine Degradation

The growth substrate can competitively inhibit micropollutant degradation because of enzymes and/or other cofactors (Fernandez-Fontaina et al., 2014) located near

U

the producing microorganism or secreted into the bulk liquid (Batstone et al., 2002). This

N

idea supports the results of the three SMZ degradation hypotheses, given that the

A

sensitivity and collinearity analyses showed that the SMZ degradation half-saturation

M

constant can be zero. Thus, the SMZ degradation is fully dependent on the OLR, from a growth and substrate consumption perspective. Table 6 shows the adjusted SMZ

ED

degradation-associated parameters, and Table 7 shows the results of the identifiability analyses. Higher OLRs corresponded to more efficient overall SMZ degradation because

PT

the hypotheses are linearly correlated with their respective substrate concentrations,

CC E

despite the decrease in the degradation/enzyme production yields. The three hypotheses showed almost the same behavior regardless of the amplitudes of the yield coefficient functions. The transformation capacity coefficients for the cometabolic hypothesis

A

indicated that approximately 0.1% of the consumed mass of the substrate was related to SMZ decomposition. In addition, during the hydrolysis stage, the three hypotheses were not sensitive to the OLR shocks, as shown in Table 7. This lack of sensitivity was similar for the SMZ shock loads, what indicates that the SMZ degradation-associated bacteria are as robust to

influent variations as the hydrolysis-related bacteria, in contrast to the acid-consuming group. This result is interesting, because the acid-consuming cometabolism was the bestfitted hypothesis, as shown in Table 4. With respect to sudden decreases in HRT, the yield sensitivity was 90% higher during the hydrolysis stage and 300% higher than the acidconsumption stage. However, positive HRT shocks were not influential for any

IP T

hypothesis. Furthermore, considering Equations (24) to (26), the optimum necessary number of days for the effluent to react to the influent variation was approximately four

SC R

days for CODt,f and SMZ and fifteen days for HRT. These are important observations, because recirculation could help avoid OLR shocks that might impact SMZ degradation. However, recirculation would require care because HRT changes also impact the process

N

U

performance, as shown during the hydrolysis stage in Figure 3.

According to Table 7, the three SMZ outputs were sensitive to νz , fνz and Yz1,2 .

A

i

M

However, evaluating these parameters in the collinearity analyses resulted in indices of 45.4, 41.7 and 57.9, respectively, for the hydrolysis and acid-consuming cometabolic

ED

hypotheses and the enzymatic degradation hypothesis. Removing any of the three parameters would cause the collinearity index to drop below 20. The first parameters

PT

reflect the microbial community adaptation rate; for most of the experiments, they would

CC E

not affect the SMZ degradation kinetics. However, if these parameters were not associated with the negative adaptation rate fraction, the decrease between E4 and E9 would not be as smooth as observed, and the degradation kinetics of those experiments

A

would be nearly unidentifiable. The SMZ transformation capacity defines the maximum degradation rate; without it, the degradation kinetics would not be identifiable. Thus, although the collinearity analysis suggested that one of these parameters should be removed, they were kept in the modeling.

Figure 6 shows the transformation capacity and the yield coefficient evolutions, and Figure 7 shows the experimental influent versus the experimental and simulated effluents. As observed in Figure 6 and Table 7, the coefficient adaptation rates were four times faster than those with substrate consumption, which is inconsistent with the degradation hypothesis. These results suggest that non-observed dynamics could be

IP T

influencing the degradation performance. The decreases in the coefficients were nearly

17 times slower than the increases caused by the OLR. This result can be explained by

SC R

the specialization of the involved bacteria group during the low OLR, which slowly decreased with increasing OLR; the same behavior was genetically observed by (Ghasimi et al., 2015). During the simulations, the adaptation rates may have been slightly lower;

U

however, these results would negatively affect the simulated effluent responses to the

N

influent variations when evaluating the profiles shown in Figure 7. During the parameter

A

estimation step, several adjustments could be performed to match the Figure 6 data.

M

However, the data in Figures 7 and 8 should be used as a reference; otherwise, the large

ED

number of degrees of freedom would result in overparameterization. Regarding the HRT effects, Figure 6 shows that a lower HRT corresponds to lower yield coefficients.

PT

As observed in Figure 7, for most of the experimental phases, the estimated

CC E

effluent SMZ was within the experimental confidence intervals. However, there were some large deviations between the estimated and measured effluent data. These deviations may stem from non-observable parts of the process, given that it was split into

A

two stages. Furthermore, given the available experimental data, increasing the number of parameters to better describe this deviation would result in parameter unidentifiability. A notable event occurred between days 280 and 350 during the organic load shock of experiment 7. The SMZ degradation hypotheses could not explain the experimental observations because the substrate concentration affected the estimates of the process

performance despite changes in the yield and transformation capacity coefficients. This finding suggests that the tested mechanisms, which were evaluated for reactions occurring near the microorganisms and for the enzymes that were secreted in the bulk liquid, were not able to represent the true process under such operational conditions. Similar behavior was observed following the sudden HRT increase between E10 and E11. For these cases,

IP T

more complex models may be required to describe the larger number of variables, as in the ADM1 (Blumensaat & Keller, 2005).

SC R

Comparing the measured SMZ at the sampling ports with the simulated hypothesis

data clarifies the strengths and weakness of each modeling structure. For experiments 2

U

to 6 (Table 1), both cometabolic degradation hypotheses represent the experimental data

N

better than the enzymatic approach, as observed in Figure 8 (a) to (e). However, after the

A

shock load of experiment 7, the enzymatic approach more accurately reflected the

M

observations (Table 4). The main advantage of the enzymatic degradation hypothesis appeared in experiments 9 and 10, in which the possibility of reactions in the bulk liquid

ED

was supported by the higher levels of SMZ along the reactor. This scenario could occur if the enzymes were washed out of the reactor because of high liquid superficial velocities.

CC E

observations.

PT

Overall, among the hypotheses, the acid-consumption hypothesis best explained the

A

5. Conclusions

In this study, a mathematical model was developed to represent the total COD

consumption, acid formation and consumption and three hypotheses for SMZ degradation. Initially, the long-term modeling considered the effects of every measured influent on the modeling. Sensitivity and collinearity analyses were conducted to avoid overparameterization, and several parameters were removed, most of them related to the

hydrolysis stage. The effects of variations in the OLR, HRT, and SMZ were quantified over the long term. The results showed that a sudden variation of 8 µg of SMZ had a greater impact on the process than did 1000 mg of filtered COD. Of the SMZ degradation hypotheses, the acid-consumption cometabolism hypothesis more accurately explained the results than did the antimicrobial removal hypothesis; however, neither was able to

IP T

accurately reflect the impacts of HRT variations on the removal process. Thus, the

enzymatic hypothesis more accurately represented the SMZ removal, which was more

SC R

pronounced after the organic shock load of the seventh experiment. This performance was especially pronounced for the lower HRT experiments, for which the enzymatic

U

hypothesis showed significantly greater accuracy compared to the other hypotheses.

N

Acknowledgments

A

The authors gratefully acknowledge the support provided for this study by the São Paulo

M

Research Foundation (FAPESP; grants 2016/15003-3, 2015/06246-7, and 2012/18942-

A

CC E

PT

ED

0).

References

A

CC E

PT

ED

M

A

N

U

SC R

IP T

Amorim, A.K.B., Zaiat, M., Foresti, E. 2005. Performance and stability of an anaerobic fixed bed reactor subjected to progressive increasing concentrations of influent organic matter and organic shock loads. Journal of Environmental Management, 76(4), 319-325. Baran, W., Adamek, E., Ziemianska, J., Sobczak, A. 2011. Effects of the presence of sulfonamides in the environment and their influence on human health. Journal of Hazardous Materials, 196, 1-15. Batstone, D.J., Keller, J., Angelidaki, I., Kalyuzhny, S.V., Pavlostasthis, S.G., Rozzi, A., Sanders, W.T.M., Siegrist, H., Vavilin, V.A. 2002. Anaerobic digestion model No. 1 (ADM1). IWA Publishing, London, UK. Blumensaat, F., Keller, J. 2005. Modelling of two-stage anaerobic digestion using the IWA Anaerobic Digestion Model No. 1 (ADM1). Water Research, 39(1), 171183. Carballa, M., Omil, F., Ternes, T., Lema, J.M. 2007. Fate of pharmaceutical and personal care products (PPCPs) during anaerobic digestion of sewage sludge. Water Research, 41(10), 2139-2150. Contois, D.E. 1959. Kinetics of bacterial growth - relationship between population density and specific growth rate of continuous cultures. Journal of General Microbiology, 21(1), 40-50. Criddle, C.S. 1993. The kinetics of cometabolism. Biotechnology and Bioengineering, 41(11), 1048-1056. Dalsenter, F., Viccini, G., Barga, M., Mitchell, D., Krieger, N. 2005. A mathematical model describing the effect of temperature variations on the kinetics of microbial growth in solid-state culture. Process Biochemistry, 40(2), 801-807. de Nardi, I.R., Zaiat, M., Foresti, E. 1999. Influence of the tracer characteristics on hydrodynamic models of packed-bed bioreactors. Bioprocess Engineering, 21(5), 469-476. Donoso-Bravo, A., Mailier, J., Martin, C., Rodriguez, J., Aceves-Lara, C.A., Vande Wouwer, A. 2011. Model selection, identification and validation in anaerobic digestion: A review. Water Research, 45(17), 5347-5364. Fernandez-Fontaina, E., Carballa, M., Omil, F., Lema, J.M. 2014. Modelling cometabolic biotransformation of organic micropollutants in nitrifying reactors. Water Research, 65, 371-383. Gabor, A., Villaverde, A.F., Banga, J.R. 2017. Parameter identifiability analysis and visualization in large-scale kinetic models of biosystems. Bmc Systems Biology, 11. Ghasimi, D.S.M., Tao, Y., de Kreuk, M., Zandvoort, M.H., van Lier, J.B. 2015. Microbial population dynamics during long-term sludge adaptation of thermophilic and mesophilic sequencing batch digesters treating sewage fine sieved fraction at varying organic loading rates. Biotechnology for Biofuels, 8. Hierholtzer, A., Akunna, J.C. 2014. Modelling start-up performance of anaerobic digestion of saline-rich macro-algae. Water Science and Technology, 69(10), 2059-2065. Horn, H., Reiff, H., Morgenroth, E. 2003. Simulation of growth and detachment in biofilm systems under defined hydrodynamic conditions. Biotechnology and Bioengineering, 81(5), 607-617. Hruska, K., Franek, M. 2012. Sulfonamides in the environment: a review and a case report. Veterinarni Medicina, 57(1), 1-35.

A

CC E

PT

ED

M

A

N

U

SC R

IP T

Kumar, S., Katiyar, N., Ingle, P., Negi, S. 2011. Use of evolutionary operation (EVOP) factorial design technique to develop a bioprocess using grease waste as a substrate for lipase production. Bioresource Technology, 102(7), 4909-4912. Larcher, S., Yargeau, V. 2012. Biodegradation of sulfamethoxazole: current knowledge and perspectives. Applied Microbiology and Biotechnology, 96(2), 309-318. Leitao, R.C., van Haandel, A.C., Zeeman, G., Lettinga, G. 2006. The effects of operational and environmental variations on anaerobic wastewater treatment systems: A review. Bioresource Technology, 97(9), 1105-1118. Michael, I., Rizzo, L., McArdell, C.S., Manaia, C.M., Merlin, C., Schwartz, T., Dagot, C., Fatta-Kassinos, D. 2013. Urban wastewater treatment plants as hotspots for the release of antibiotics in the environment: A review. Water Research, 47(3), 957-995. Mohring, S.A.I., Strzysch, I., Fernandes, M.R., Kiffmeyer, T.K., Tuerk, J., Hamscher, G. 2009. Degradation and Elimination of Various Sulfonamides during Anaerobic Fermentation: A Promising Step on the Way to Sustainable Pharmacy? Environmental Science & Technology, 43(7), 2569-2574. Oliveira, G.H.D. 2016. Remoção de sulfametazina em reatores anaeróbios tratando água residuária de suinocultura. in: Escola de Engenharia de São Carlos, Universidade de São Paulo, pp. 193. Oliveira, G.H.D., Santos-Neto, A.J., Zaiat, M. 2016. Evaluation of sulfamethazine sorption and biodegradation by anaerobic granular sludge using batch experiments. Bioprocess and Biosystems Engineering, 39(1), 115-124. Oliveira, G.H.D., Santos-Neto, A.J., Zaiat, M. 2017. Removal of the veterinary antimicrobial sulfamethazine in a horizontal-flow anaerobic immobilized biomass (HAIB) reactor subjected to step changes in the applied organic loading rate. Journal of Environmental Management, 204, 674-683. Onesios, K.M., Yu, J.T., Bouwer, E.J. 2009. Biodegradation and removal of pharmaceuticals and personal care products in treatment systems: a review. Biodegradation, 20(4), 441-466. Perez, S., Eichhorn, P., Aga, D.S. 2005. Evaluating the biodegradability of sulfamethazine, sulfamethoxazole, sulfathiazole, and trimethoprim at different stages of sewage treatment. Environmental Toxicology and Chemistry, 24(6), 1361-1367. Pomies, M., Choubert, J.M., Wisniewski, C., Coquery, M. 2013. Modelling of micropollutant removal in biological wastewater treatments: A review. Science of the Total Environment, 443, 733-748. Ramirez, I., Mottet, A., Carrere, H., Deleris, S., Vedrenne, F., Steyer, J.P. 2009. Modified ADM1 disintegration/hydrolysis structures for modeling batch thermophilic anaerobic digestion of thermally pretreated waste activated sludge. Water Research, 43(14), 3479-3492. Raue, A., Kreutz, C., Maiwald, T., Bachmann, J., Schilling, M., Klingmueller, U., Timmer, J. 2009. Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics, 25(15), 1923-1929. Saravanan, V., Sreekrishnan, T.R. 2006. Modelling anaerobic biofilm reactors - A review. Journal of Environmental Management, 81(1), 1-18. Sarti, A., Vieira, L.G.T., Foresti, E., Zaiat, M. 2001. Influence of the liquid-phase mass transfer on the performance of a packed-bed bioreactor for wastewater treatment. Bioresource Technology, 78(3), 231-238.

A

CC E

PT

ED

M

A

N

U

SC R

IP T

Shelver, W.L., Hakk, H., Larsen, G.L., DeSutter, T.M., Casey, F.X.M. 2010. Development of an ultra-high-pressure liquid chromatography-tandem mass spectrometry multi-residue sulfonamide method and its application to water, manure slurry, and soils from swine rearing facilities. Journal of Chromatography A, 1217(8), 1273-1282. Singh, B., Kumar, R., Ahuja, N. 2005. Optimizing drug delivery systems using systematic "Design of experiments." Part I: Fundamental aspects. Critical Reviews in Therapeutic Drug Carrier Systems, 22(1), 27-105. Vera, D.R., Scheibe, P.O., Krohn, K.A., Trudeau, W.L., Stadalnik, R.C. 1992. GOODNESS-OF-FIT AND LOCAL IDENTIFIABILITY OF A RECEPTORBINDING RADIOPHARMACOKINETIC SYSTEM. Ieee Transactions on Biomedical Engineering, 39(4), 356-367. Villaverde, A.F., Banga, J.R. 2014. Reverse engineering and identification in systems biology: strategies, perspectives and challenges. Journal of the Royal Society Interface, 11(91). Wang, N., Guo, X., Xu, J., Kong, X., Gao, S., Shan, Z. 2014. Pollution characteristics and environmental risk assessment of typical veterinary antibiotics in livestock farms in Southeastern China. Journal of Environmental Science and Health Part BPesticides Food Contaminants and Agricultural Wastes, 49(7), 468-479. Wang, Y.J., Witarsa, F. 2016. Application of Contois, Tessier, and first-order kinetics for modeling and simulation of a composting decomposition process. Bioresource Technology, 220, 384-393. Yin, X., Qiang, Z., Ben, W., Pan, X., Nie, Y. 2014. Biodegradation of Sulfamethazine by Activated Sludge: Lab-Scale Study. Journal of Environmental Engineering, 140(7).

Set the initial values for the simulations in experiment 7 equal to the end values and consider the biomass to be fully acclimatized to the OLR, e.g., F1,2 = 1.0.

Adjust the k1 , k 2 , Yh ,Ya , K d and Ds values to adjust the simulated CODf, VFA and biomass values to the measured values; otherwise, adjust only the k1 and k 2 values.

IP T

Hold the Yh ,Ya , K d and Ds values constant, and adjust k1 and k 2 to approximate the simulated CODf and VFA values for experiments 1 to 7.

Are the optimum ajustment and Chisquare bellow critical value?

No

Yes

SC R

Adjust 𝑘1 , 𝑘 2 , 𝐾𝑑 and 𝐷𝑠 to adjust the simulated CODf and VFA values for experiments 8 to 11

Is the adjustment value better than previous optimum value?

N

No

A

Hold the Yh ,Ya , K d and Ds values.

U

Use distinct Yh ,Ya , K d and Ds values to adjust the biomass, CODf and VFA in experiment 11.

Yes

M

Select SMZ degradation

A

CC E

PT

ED

Given the mechanism, hold the k1 , k 2 , Yh ,Ya , K d, and Ds values constant and adjust the YE and K De or Ys values for experiments from 1 to 11.

Figure 1: Fitting procedure flowchart.

Given the static parameters, adjust Equations (18) to (23) to the calculated k1 , k 2 , YE and Ys values subject to the OLR and HRT.

Is the chi-square values bellow critical value? Yes End

No

IP T SC R U N A

M

Figure 2: Influent and effluent COD data, respectively. In (a), the thick black line and thin gray line are the experimental CODt and experimental CODf influents, respectively. In

ED

(b), the thin blue line is the measured soluble effluent, and the thick blue line is the

PT

simulated soluble effluent. The thick and the thin dashed lines represent the confidence

A

CC E

interval and the mean of the experimental effluent COD, respectively.

IP T SC R

U

Figure 3: Comparison between the substrate consumption rate parameters K1 and K2 for

A

CC E

PT

ED

M

A

N

both stages under OLR and HRT variations with their simulated behaviors K1m and K2m.

Figure 4: Experimental and simulated data for VFAs and CODf at the sampling ports. The red ‘x’ and yellow ‘o’ marks indicate the experimental CODf and VFA observations, respectively. The blue lines represent the simulated CODf consumption, the black lines represent the VFA dynamics, and the dark green lines represent the effluent CODf confidence intervals, normalized and leveled for each sampling port. The subplots show

IP T

the following experiments: (a) E2, (b) E3, (c) E4, (d) E5, (e) E6, (f) E7, (g) E6, (h) E9,

M

A

N

U

SC R

(i) E10 and (d) E11.

ED

Figure 5: Comparison between the measured and simulated biomass for (a) E7 and (b)

A

CC E

PT

E11.

IP T SC R U

N

Figure 6: Sulfamethazine degradation coefficient profiles. The ‘x’ marks and dashed

A

black lines show the adjusted parameters from Table 4; yellow, purple and green lines

M

indicate the enzymatic production yield and the hydrolysis-related and acid-consumption-

A

CC E

PT

ED

related cometabolic capacities in (a), (b) and (c), respectively.

Figure 7: Comparison between the influent and effluent SMZ concentrations during the operational phases. The blue and red lines indicate the experimental influent and effluent concentrations, respectively, and the yellow, purple and green lines represent the simulated data for the enzymatic, hydrolysis-related and acid-consumption-related

cometabolic degradations, respectively. The thick and thin dashed lines represent the

ED

M

A

N

U

SC R

IP T

confidence intervals and the mean for each operational phase, respectively.

Figure 8: Experimental and simulated SMZ data at the sampling ports. Black ‘x’ marks

PT

represent the experimental observations; the blue, red and dashed yellow lines represent

CC E

the enzymatic, acid-consumption-related and hydrolysis-related degradation simulated kinetics, respectively. Subplots show experiments (a) E2, (b) E3, (c) E4, (d) E5, (e) E6,

A

(f) E7, (g) E6, (h) E9, (i) E10 and (j) E11.

Table 1: Operating conditions of each experiment (Oliveira et al., 2017) Operational Phase

Duration (days)

HRT (h)

Average 𝐂𝐎𝐃𝐭

Average 𝐂𝐎𝐃𝐟

(𝐤𝐠 𝐎𝟐 . 𝐦−𝟑 . 𝐝−𝟏 )

(𝐤𝐠 𝐎𝟐 . 𝐦−𝟑 . 𝐝−𝟏 )

E1

35

24

3.0±0.2

2.1±0.1

0.0

E2

42

24

3.0±0.2

2.1±0.1

8.8±0.4

E3

66

24

1.8±0.1

1.1±0.2

8.7±0.9

E4

29

24

0.9±0.1

0.6±0.1

7.4±0.7

E5

33

24

2.3±0.4

1.5±0.3

9.0±0.8

E6

40

24

2.9±0.4

2.1±0.2

8.8±0.4

E7

25

24

6.2±0.8

3.8±0.7

9.6±1.0

E8

19

24

3.1±0.5

2.1±0.4

10.3±0.4

E9

26

16

5.0±0.6

3.1±0.1

13.5±0.6

E10

18

8

9.8±1.0

6.2±0.5

28.2±0.8

E11

78

24

3.1±0.8

1.4±0.2

9.6±0.6

PT

ED

M

A

Biomass Growth and Substrate Consumption CODf and Analysis SMZC SMZE VFA NP 6 7 8 NS 12 18 18 df 5 10 9 CV (95%) 11.07 18.3 16.9

CC E

IP T

SC R

U

N

Table 2: Chi-square critical values for goodness-of-fit.

A

Average SMZ (𝐠 𝐒 𝐦𝐳 . 𝐦−𝟑 𝐝−𝟏 )

Table 3: Optimal adjustments of fixed and varying parameters. Definition

Estimated Value

Yh

Hydrolysis-related biomass growth yield coefficient

7.0× 10−6

Ya

Acid-consumptionrelated biomass growth yield coefficient

K Sh

Hydrolytic biomass concentration halfsaturation constant

K Sa

Acid-consumption biomass concentration halfsaturation constant

KD

Endogenous biomass decay rate constant

DS

Sludge time

F1

Hydrolytic biomass adaptation factor

1.0*

F2

Acid-consumption biomass adaptation factor

1.0*

in

E2

20 00

IP T

mg CODf . g F . g −1 SSV

mg CODf . g F . g −1 SSV

SC R

830

h−1 h−1

N

3.5 × 10−6

U

0

A

-

E9

E10

E11

20 00

20 00

2000

1200

36. 4

36. 4

47. 0

60. 6

30. 3

38. 6

30. 2

41.4

27.0

77 0

440

92 4

15 40

40. 7

44. 4

82. 0

117.1

230

7.9

15.8

6.3

PT

CC E

E8

40 00

53. 0

χ2VA

E7

20 00

−1 mg CODf . g F L−1 g −1 SSV cm 3 × 10

χ2CODf

E6

15 00

K 2m

𝐑𝟐𝐕𝐀

E5

540

40. 9

𝐑𝟐𝐂𝐎𝐃𝐟

E4

10 00

−1 mg CODf . g F L−1 g −1 SSV cm33. 3 0 × 10

x10−6 h−1

E3

20 00

K 1m

A

−1 g SSV . L. g −1 F . mg CODf

830

M

E1

(mg O2 . L−1 . d−1 )

DS

−1 g SSV . L. g −1 F . mg CODf

Experimental Conditions

Parameters COD f

Unit

2.2 × 10−6

retention

ED

Fixed Parameters

Parameters

3.50 0.9 9

0.9 9

0.9 9

0.9 5

0.9 9

0.9 9

0.9 9

0.9 9

0.9 9

0.99

0.99

-

0.9 7

N C

NC

N C

N C

0.9 9

0.9 6

0.9 4

0.89***

0.06* *

4.6 6

2.4 3

2.4 2

27. 78

0.5 9

1.6 1

1.3 6

1.2 6

0.7 1

3.58

0.23

-

0.4 8

0.0 8

0.0 2

0.0 1

0.0 1

1.3 8

3.3 6

1.5 5

2.16***

0.11

* Assumed ** This correlation is not valid, because of problems with the R2 calculation method; the experimental observations that are close to but not equal to zero, as shown in Figure 3.

*** For the R2 and for chi-square calculations, only valid measurements were considered; the first sampling device had missing data. In this case, the 𝛘𝟐 critical value was 9.48.

Table 4: Sulfamethazine degradation-adjusted parameters

R2Sh R2Sa R2E 𝜒𝑆2ℎ 𝜒𝑆2𝑎

A

CC E

PT

𝜒𝐸2

E5

E6

E7

E8

E9

E10

E11

320

800

1380

1075

695

300

473

250

140

680

345

800

1400

1050

700

365

510

245

140

690

158

365

650

515

335

140

210

92

65

300

0.90

0.88

0.98

0.99

0.90

0.88

0.98

0.89

0.87

1.64

1.82

1.61 1.52

∑ 𝜒2

SC R

IP T

E4

0.99

0.80

0.77

0.32

0.19

0.91

N

YE

E3

0.99

0.99

0.94

0.98

0.61

0.43

0.93

A

YSa

Hydrolysisrelated cometabolic transformation capacity Acidconsuming μg Smz cometabolic mg CODf transformation capacity Hydrolysisrelated U enzymatic mg CODf production yield constant Hydrolysis-related cometabolism R2 Acid-consuming cometabolism R2 Apparent enzymatic activity degradation R2 Hydrolysis-related cometabolism χ2 Acid-consuming cometabolism χ 2 Apparent enzymatic activity degradation χ 2 μg Smz mg CODf

E2

0.93

0.97

0.92

0.96

0.98

0.76

0.73

0.98

0.65

0.61

1.25

3.01

5.36

5.68

10.8

0.69

31.5

1.80

0.61

0.64

1.19

0.79

1.00

3.01

6.25

0.58

17.5

1.94

3.34

2.08

6.99

0.89

1.02

1.87

2.95

0.27

22.9

M

YSh

Unit

ED

Parameter

U

Experimental Conditions Definition

Table 5: Parameters associated with long-term OLR and HRT variations

K Sz 1

Hydrolysis Stage

K v1 fv1 n m K 2m K s2 ν2

CC E

K e2 p

Acid-Consumption Stage

f2

A

K Sz 2 K v2 fv2

n m

630

mg CODf . g F . L−1 g −1 SSV mg CODf . L−1

630

1.95 × 10−3

1.70 × 10−3

624.1 × 103

0

h−1

IP T

69.7 × 103

Unit

mg 2CODt . L−2 . h

2

-

-

0.03

0

4.20

∞

0.79

SC R

f1

66.7 × 103

U

p

After S/C Analyses

-

μg 2SMZ . L−2 h

0.79

L

0.38

−

4

-

d

11

11

d

N

K e1

Estimated Value

0.75

M

ν1

ED

K s1

Maximum hydrolysis/acid formation capacity constant Half saturation for hydrolytic/acid formation bacteria constant Hydrolytic/acid formation bacteria adaptation rate constant Hydrolytic/acid formation bacteria OLR shock sensitivity constant Hydrolytic/acid formation bacteria OLR shock sensitivity constant Hydrolytic/acid formation bacteria Negative OLR sensitivity fraction constant Hydrolytic/acid formation bacteria SMZ shock sensitivity constant Hydrolytic/acid formation bacteria negative HRT shock sensitivity constant Hydrolytic/acid formation bacteria positive HRT shock sensitivity fraction constant Process response time to OLR and SMZ variations Process response time to HRT variations Maximum acid-consumption capacity constant Half saturation for acidconsumption bacteria constant Acid-consuming bacteria adaptation rate constant Acid-consuming bacteria OLR shock sensitivity constant Acid-consuming bacteria OLR shock sensitivity constant Acid-consuming bacteria negative OLR sensitivity fraction constant Acid-consuming bacteria SMZ shock sensitivity constant Acid-consuming bacteria negative HRT shock sensitivity constant Acid-consuming bacteria positive HRT shock sensitivity constant Process response time to OLR and SMZ variations Process response time to HRT variations

PT

K 1m

Definition

A

Parameter

3.64 × 106 960

3.64 × 106

mg CODf . g F . L−1 g −1 SSV mg CODf . L−1

960

2.0 × 10−3

1.9 × 10−3

h−1

133.6 × 103

130.4 × 103

mg 2CODt . L−2 . h

2

2

-

0.03

0.047

-

1.43

1.15

μg 2SMZ . L−2 h

1.28

1.28

L

0.46

1

−

4

4

d

11

15

d

Table 6: Results of sensitivity and collinearity analysis of parameters Hydrolysis and acid-consumption related parameters Ya

K 1M

K 2M

K s1

K s2

K s𝑧1

K s𝑧2

𝐾𝑠𝑋

DS

ν1

ν2

K e1

CODf

S/C

S/C

S/C

S/C

S/C

S/C

S

S/C

S/C

S/C

S/C

S/C

S

VFA

S/C

S/C

S/C

S/C

S/C

S/C

S/C

S/C

S/C

S/C

S/C

S/C

S/C

X

S

S

S

S

S

S

S

S

S

S/C

S

S

SmzH

S/C

C

C

S/C

C

C

C

S/C

C

S/C

C

S/C

C

SmzA

C

C

C

C

C

C

C

S/C

C

C

C

S/C

C

SmzE

C

C

C

C

C

C

C

S/C

C

S/C

C

S/C

C

f1

K e2

f2

K v1

fv1

K v2

fv2

n

m

S

S/C

S/C

S/C

C

S/C

C

S/C

S/C

S/C

S/C

S/C

S/C

C

S/C

C

S/C

S/C

S

S/C

S

S/C

S/C

S/C

S/C

S/C

C

C

C

S/C

S/C

C

C

S/C

S/C

C

C

S

C

C

U N A M ED PT CC E A

C C

C C

S

C

C

C

C

S

C

IP T

Yh

SC R

Output

Table 7: SMZ degradation-associated parameters with OLR and HRT variations

fνz1

fνz2

Cometabolic SMZ degradation

K 𝑠z1

M

ED

K sz

CC E

n

PT

K v𝑧 fv

SMZ

m

YE

Bulk enzymatic degradation

A

K zE νz

fν z

E

0.78

μg SMZ . mg −1 COD

50

0

8.0 × 10−3

0.02

0.02

f

mg COD . L−1 f

8.0 × 10−3

0.069

0.0645

h−1

−

−

-

mg 2CODt . L−2 . h

523.9 × 103

-

mg 2CODt . L−2 . h

8.29

-

mg 2SMZ . L−2 h

0.51

0.41

L

1.59

1

−

4

4

d

11

15

d

0.39

0.36

50

0

523.9 × 103

A

K sz2

0.89

IP T

νz

Unit

SC R

K z1,2

SMZ transformation capacity for hydrolysis cometabolism SMZ transformation capacity half-saturation constant associated with OLR SMZ degradation-associated bacteria adaptation rate constant SMZ degradation-associated bacteria negative adaptation rate constant fraction for hydrolysis cometabolism SMZ degradation-associated bacteria negative adaptation rate constant fraction for acids consuming cometabolism Hydrolysis relate SMZ transformation capacity sensitivity constant to OLR shock loads Acids consuming relate SMZ transformation capacity sensitivity constant to OLR shock loads SMZ transformation capacity sensitivity constant to SMZ shock loads SMZ transformation capacity sensitivity constant to positive HRT variations SMZ transformation capacity sensitivity constant to negative HRT variations Process response time to OLR and SMZ variations Process response time to HRT variations Maximum enzymatic production yield Enzymatic production yield half-saturation constant Enzymatic production associated bacteria adaptation rate constant fraction Enzymatic production associated bacteria negative adaptation rate constant fraction for enzymatic degradation

After S/C Analyses

U

Yz1,2

Estimated Value

N

Parameter Definition

8.0 × 10−3

0.02

8.0 × 10−3

0.069

μg SMZ . mg −1 CODf mg CODf . L−1 h−1

−

K sz K vz fv K De n m

523.9 × 103

-

mg 2CODt . L−2 . h

8.29

-

μg 2SMZ . L−2 h

0.51

0.41

L

1.59

0

−

0.40

0.40

4

4

11

ν z1 S/C* C C

ν z2 C S/C* C

ν ze C C S/C*

fvz1 S/C* C C

fvz2 C S/C* C

d

K 𝑧1,2,𝐸

C C S/C*

K s𝑧1 S/C* C C

K s𝑧2

C S/C* C

K s𝑧 𝐸 C C S/C*

fz1,2,E

K sz

K v𝑧 S S S

CC E

PT

ED

M

*By some measures, these parameters were considered influential to the output.

A

d

U

fvzE

A

SmzH SmzA SmzE

YE C C S/C

N

Sulfamethazine degradation-related parameters Yz2 C S/C C

h−1

15

Table 8: Sensitivity and collinearity analysis results

Yz1 S/C C C

IP T

E

SC R

Enzymatic production yield sensitivity constant to OLR shock Enzymatic production yield sensitivity constant to SMZ shock loads Enzymatic production yield sensitivity constant to positive HRT variations Enzymatic production yield sensitivity constant to negative HRT variations Apparent enzymatic activity degradation rate Process response time to OLR and SMZ variations Process response time to HRT variations

K sz

fv

K De C C S/C

n

m

S/C S/C S/C

S/C S/C S/C