Modelling Emerging Pollutants in Wastewater Treatment: A Case Study using the Pharmaceutical 17 α − ethinylestradiol

Mathematical modelling can play a key role in understanding as well as quantifying uncertainties surrounding the presence and fate of emerging pollutants in wastewater treatment processes (WWTPs). This paper presents for the ﬁrst time a simpliﬁed emerging pollutant pathway in the WWTP that incorporates two potential pathways to sequestration. It develops deterministic and stochastic ordinary diﬀerential equations to gain insight into the fate and behaviour of a case study pharmaceutical, with particular focus on sorption to the solid phase, as well as the nature of the experimentally measured solid parent compound. Statistical estimation and inferential pro-cedures are developed and via a proof-of-concept examination, the study explores the transformation pathways of the bioactive chemicals (BACs) in the bioreactor, which is the heart of the WWTP. With a focus on the case study pharmaceutical 17 α − ethinlyestradiol (EE2), the simulation results show good agreement with the EE2 data. In addition, the results suggest that the experimentally measured solid EE2-parent concentration is very similar to the model-based sequestered EE2-parent concentration.


Introduction
Maintaining water quality is a major global challenge due to the increasing use of pharmaceutical drugs and related bioactive chemicals (BACs).
Since all water is reused, the urban water cycle requires protection via treatment to ensure that water is of appropriate quality. However, wastewater treatment processes (WWTPs) were never designed to remove emerging pollutants (such as the BACs) from wastewater (Margot et al., 2015;Salgado et al., 2012). Therefore, in order to improve BAC removal it is important to understand the transformation pathways of BACs in the WWTP.
The presence of BACs (e.g., pharmaceutical, steroids) observed in wastewater and the aquatic environment even at low concentrations (ngL −1 ) can cause adverse effects in aquatic organisms (Verlicchi et al., 2012). Hence, increasingly stringent legislation such as the EU Water Framework Directive (WFD) 2015/495 requires those responsible for wastewater treatment to monitor the presence of these pollutants with the future aim of improving the quality of the wastewater effluent released to surface waters (Carvalho et al., 2015;Barbosa et al., 2016). This has resulted in a significant number of investigations into the removal of BACs in the WWTP. Recent studies include mathematical models that take into consideration the influence of free (or parent) compounds as well as retransformable products such as conjugated compounds (Xu et al., 2016), where conjugated compounds are Phase II metabolites formed as a result of drug transformation processes (Coleman, 2010).
Sorption to activated sludge plays a role in the fate of BACs in WWTPs and is dependent on the physiochemical properties of the BAC. Sorption can be a precursor to biotransformation/biodegradation. Several experimental sorption studies have found the occurrence of sequestration, where pollutants associated with the solid phase (e.g. mixed liquor suspended solids (MLSS), sludge or sediment) is not readily reversible (Yi and Harper, 2007;Barret et al., 2011;Sittig et al., 2012). However, there is limited mathematical modelling insight into this phenomenon and no studies to the best of the authors knowledge show whether the BACs concentration on the solid phase is on the surface of the activated sludge (adsorption), within (absorption) or a combination of both (Gomes et al., 2009;Banihashemi and Droste, 2014).
Distinguishing between adsorption and absorption is vital since the absorbed BACs are not readily available for microbial biodegradation (Banihashemi and Droste, 2014). That is, the ratio of BACs adsorbed to absorbed has implications in terms of the distribution from liquid phase to solid phase and vice versa (Yi and Harper, 2007;Barret et al., 2011). This has not been the focus of many investigations but as Barret et al. (2011) suggested, the quantification of sorption of BACs with activated sludge will aid the prediction of the fate of BACs and the associated risk to the environment.
The majority of mathematical modelling research carried out in the field of WWTP (lab-, pilot-or full-scale) employs deterministic ODE modelling based on the assumption that all the biological, chemical and physical laws governing the system can be understood (Gernaey et al., 2004). For example, mathematical models describing the fate of both traditional and emerging pollutants in the WWTP have been placed under the general framework of the activated sludge model for Xenobiotics (ASM-X) (Plósz et al., 2010;Snip et al., 2014;Polesel et al., 2015). Nevertheless, modelling of WWTP is adversely influenced by imperfect knowledge of biochemical processes, operational parameters, environmental factors and input loadings. These variabilities have both temporal and spatial effects in the removal of emerging pollutants from treated wastewater. Liu et al. (2015) pointed out that while lab-scale batch experiments have a role, they do not often reflect fullscale wastewater treatment or integrated wastewater treatment observations because wastewater is stochastic in nature. Stochastic models allow the estimation of model parameters and the understanding of the uncertainties in pollutant concentrations.
EE2 is a pharmaceutical and synthetic estrogen used in birth control pills (Laurenson et al., 2014). A review conducted by Ting and Praveena (2017) reveals that the excretion rate for women who used such medication is 35µgday −1 and that 80% of EE2 taken into the body is excreted as unmetabolized conjugates with 30% from faeces and 22-50% from urine. Due to the potential adverse endocrine disrupting effects associated with even small (ngL −1 ) concentrations in the aquatic environment, it is reported in the Watch List of the EU WFD 2015/495 (Carvalho et al., 2015;Barbosa et al., 2016). Thus, EE2 is considered as the case study BAC in this paper and the pathway of EE2 in WWTPs proposed by Joss et al. (2006), which is widely applied to model the fate of most emerging pollutants, is extended by introducing the concept of sequestration to investigate the relationship and nature of BACs and the solid phase. The impetus for this research is the fact that most work on sorption used an experimental (Yi et al., 2011) rather than a mathematical modelling approach. It should be noted that although Polesel et al. (2015) and Snip et al. (2014) considered the formation of sequestered products in their full-scale wastewater treatment modelling for other pharmaceuticals, they assumed that the ratio of sequestered concentration to pre-clarified influent concentration of the free parent compound is constant and as such, no kinetic rate parameter is accounted for the sequestration mechanism.
Mathematical modelling saves time and the cost of repeatedly designing new experiments by offering the opportunity to alternatively deepen knowledge of complex mechanisms via computer simulation evaluations. For this reason, a mathematical modelling approach using experimental batch data taken from an activated system (Gomes et al., 2009) is used with the aim of extracting the maximum information on the mechanisms involved in EE2 treatment by fitting a range of advanced mathematical models, including both deterministic, stochastic and combined approaches, to a dataset. This study is concerned with four forms of EE2, which are: the parent chemical in the liquid phase (free); adsorbed onto activated sludge (sorbed), absorbed into activated sludge (sequestered) and those retransformable forms of the parent chemical (e.g., excreted Phase I and II metabolites) as summarised in Table 1 of Section 3. An experimental dataset from Gomes et al. (2009) is selected for the study since it captures the changes in the concentrations of EE2 states over time. Data of such a nature are limited in the literature as more usually data involve monitoring influent and effluent only (Ting and Praveena, 2017;Yu et al., 2019).

Experimental data
The data are obtained from batch scale fate studies using activated sludge from a pilot plant following a Husmann design (Gomes et al., 2009).
A detailed description of the Husmann pilot plant used by Gomes et al. (2009) follows the OECD guidelines for testing chemicals in an activated sludge process (OECD, 2001). The dataset contains measurements of concentrations at eight different time points for three EE2 forms: the first form is the aqueous parent EE2, the second is the aqueous conjugated steroid ethinylestradiol 3 -glucuronide (EE2 -3G) and the third one is the solid parent EE2. The measurements are recorded in hours (h) at 0, 0.0167, 0.5, 1, 2, 4, 8 and 24 hours in triplicate. The mixed liquor suspended solid concentration (C MLSS ) reported is 4 gL −1 (Gomes et al., 2009). For a detailed description of the experimental setup see Gomes et al. (2009). For the purpose of mathematical modelling, Y is used to denote the measured EE2 concentrations, where Y 1 is the liquid form of parent EE2, Y 2 is the EE2-3G and Y 3 is the solid form of parent EE2.

Model description
A computational model was developed to describe the fate of BAC(s) in the unit treatment technology utilising lab-scale data. In experimental studies, the activated sludge is usually spiked with the BAC(s) of interest, which are classified by four state variables, namely, the liquid form of parent compound, C f ree , liquid form of retransformable compounds, C con , sorbed parent compounds, C sor , and sequestered parent compounds, C seq (Polesel et al., 2015). The activated sludge bioreactor is assumed to contain a constant biomass (MLSS) denoted by C MLSS = 4 gL −1 (Gomes et al., 2009) since the BAC concentrations are in the range of mgL −1 to ngL −1 , the biomass growth is negligible during the consumption of the BACs and any effects of the BAC on the actual biomass activity are assumed to be negligible (Snip et al., 2014). The following additional assumptions are made: 1. the reactor, of volume V , is well-mixed and operates on a continuously stirred tank reactor (CSTR), 2. the dissolved oxygen concentration in the reactor is constant, 3. the retransformable compounds (C con ) are not further biotransformed, but may be deconjugated back to the liquid form of parent compounds (C f ree ) (Joss et al., 2006;Gomes et al., 2009;Kumar et al., 2012;Snip et al., 2014), 4. the dissolved forms of parent compounds (C f ree ) in liquid are in a reversible reaction with the solid (sorbed) fraction (C sor ) and are further biodegraded via microbial activities as well as irreversibly absorbed, i.e., sequestered to C seq . 5. C MLSS is assumed constant. Figure 1 shows the simplified pathways diagram, which is modelled using first-order and pseudo-first order kinetics. Table 1 gives a summary of the notation, description of the model parameters and their corresponding units.  Table 1 for an explanation of parameters and concentration variables.   Figure 1 depicts the interaction between C f ree , C con , C sor and C seq , it can be seen that microbial biodegradation and sorption activities are the two main mechanisms through which the EE2-parent can be removed. Evolution of the four concentration forms of EE2 within the bioreactor can be described by four coupled nonlinear systems of ODEs given as With respect to the lab-scale batch experimental set-up, one state variable has non-zero initial concentration, C con (0) = C con 0 = 2125 ngL −1 (Gomes et al., 2009) and the remaining state variable concentrations are all zero (C f ree (0) = 0, C sor (0) = 0 and C seq (0) = 0). Note that M f ree (g mol −1 ) and M con (g mol −1 ) are respectively the molecular masses of C f ree and C con and the fraction M f ree /M con corrects for the change from one compound form to another.
The novelty of the proposed model given by Equations ( (sequestered) parent compound concentration (C seq ). However, in the proposed model, the kinetic rate parameters k ab,f ree and k ab,sor account for the absorption of both C f ree and C sor respectively. The ratio of BACs adsorbed to absorbed have implications for the actual fate and the availability to and from the liquid and solid phases (Yi and Harper, 2007). As a result, the quantification of sequestered BACs will offer a useful insight and help reduce the variability surrounding the fate of BACs in wastewater treatment.
To take into consideration the fact that experimental data are likely to have no sequestration or variation in sequestration pathways since experimental data rarely distinguishes between adsorbed and absorbed BACs.
Therefore, modifications to the proposed pathways ( Figure 1) are made leading to M0, M1 and M2 as depicted in Figure 2. The pathways (1), (2), (3) and (4) are well established pathways used in modelling the fate of BACs in wastewater treatment (Joss et al., 2006;Plósz et al., 2012b;Polesel et al., 2015;Plósz et al., 2012a) and so is included in all three models. Thus, the three different sets of pathways this study wish to compare are defined as follows: • M0, where all pathways are included, • M1, where sequestration via sorbed parent BACs is neglected, • M2, where sequestration via liquid parent BACs is neglected. The ODE bioreactor model of BACs as given by equations (3.1) -(3.4) and modifications can be written in the general framework given by the following four-dimensional system of ODEs for the underlying state variable, X t , as with X 0 = [0, C con,0 , 0, 0] T as initial state concentrations. Note that, t denote time; (3.6) are the concentrations at time t; T superscript denote matrix transposition; is the rate parameter vector; and f (X t , u t , t, θ 1 ) is a nonlinear function with f (·) describing the nonlinear relationship between the underlying state variable, X t and other input parameters, i.e., u t = C MLSS as the input parameter. For simplicity, in the method below the notation X t and θ 1 are rather used than that of Table 1.

Stochastic ODE model
The state of the biological treatment phase of the WWTP is affected by a number of factors (temperature, pH, product formation, mixing, etc.), which cannot be precisely accounted for by the ODE model (3.5) (Sathyamoorthy et al., 2014). As a result, processes such as conjugation, deconjugation, biotransformation, sorption, desorption, are affected by stochastic fluctuations in the concentrations due to (i) imprecise knowledge of rate coefficients, θ in (3.5) (e.g., dependency on temperature, pH, etc.) (ii) structural misspecification, f (·) in (3.5) due imperfect knowledge of the biological, chemical and physical processes (e.g., some compounds/products acting as inhibitors to other reactions) and (iii) fluctuations in variable concentrations due to, for example, imperfect mixing as well as varying load entering the treatment process.
To account for the uncertainties and simplifications in the previous model given by equations (3.1) -(3.4), an extension to the SDEs is considered by incorporating multiplicative system noise into the system by equations (3.1) t denotes a standard Wiener process, also known as Brownian motion (Øksendal, 2003). The system given by equations (3.8) -(3.11) can be written in a general framework by the following four-dimensional system of SDEs for the underlying state variable, X, as ( 3.12) where f (X t , u t , t, θ 1 )dt denotes the drift term with f (·); representing the deterministic part; σ(X t , u t , t, θ 1 )dW t denotes the stochastic part of the system with and θ 1 = [k dc , k d , k de , k bio , k ab,f ree , k ab,sor , σ f ree , σ con , σ sor , σ seq ] denotes the vector of model parameters in the dynamic model (3.12). All other notation remains the same as previously defined in Section 3.2.1. It is worth mentioning that the ODE model (3.1-3.4) is the SDE model (3.8-3.11) with σ f ree = 0, σ con = 0, σ sor = 0, σ seq = 0.

Model for the observations
It is often the case that in wastewater treatment not all components of the system can be observed. Hence, the observational model is vital to statistically assess the disparity between the observed data and the model prediction. The observational model relates the dataset of external observations, Y , to the internal state variables, X. As described in section 2, three forms of EE2 measurements are available for this study. However, for the measurement with respect to EE2-parent concentration in solid phase, Y 3 it cannot be differentiated as to whether Y 3 represents: H1. the sorbed concentration, C sor , H2. the sequestered concentration, C seq , H3. the sum of a fraction of the sorbed parent compound together with the sequestered parent compound, f sor C sor + C seq with 0 ≤ f sor ≤ 1.
The aim is to elucidate which of these three hypotheses (H1, H2 and H3) is most appropriate by fitting the data Y to the model (3.5) or (3.12).
However, consideration shall only be given to model (3.12), since it is a more general model as discussed earlier in 3.2.2. The observational model can be written in the general framework as: where t n with n = 1, 2, · · · , N are the discrete time points; (q = 1, 2, 3 = Q) is the number of repeated set; Y n,q = [Y 1 , Y 2 , Y 3 ] T n,q are the experimentally measured concentrations at time point t n for the q th experiment; θ 2 is the observational model parameter which will be defined in subsequent text and h(·) denotes the observational function which we define here by equation (3.15) (see Text S1 in the supporting information for further details). Note, The unknown independent measurement errors at time point t n for the q th experiment, e n,q are assumed to belong to a multivariate normal distribution with mean zero and covariance matrix given by either • or different measurement variances, The covariance structure of the measurement noise given by equation (3.16) or (3.17) is called additive noise and it is the most fundamental noise structure considered for this study due to the nature of the experimental data (i.e. the observations at different times are assumed to have a common vari-ance). Furthermore, the experimentally observed EE2 concentrations Y 1 , Y 2 and Y 3 are assumed to be uncorrelated in all triplicate experimental data because of the experimental design used.

Estimation model parameters
The aim is to fit θ = [k dc , k d , k de , k bio , k ab,f ree , k ab,sor , f sor , σ f ree , σ con , σ sor , σ seq , σ 2 0 , σ 2 1 , σ 2 2 , σ 2 3 ] in the more general model which is the stochastic ODE as described in Sections 3.2.2 and 3.2.3. In this paper the method of parameter estimation of SDE proposed by Picchini et al. (2006) is applied to the problem under study, that is, equations (3.12) and (3.18). This is a method used to carry out statistical inference in discretely observed diffusion process with partially observed state variables. From (3.18), the simulated maximum likelihood function approximation of the conditional density function for each hypothesis, ϕ = {1, 2, 3} is where, R denotes the number of Monte Carlo simulations, X (r) n denotes the r th simulated solution, X n of the equation (3.12) using the Euler-Maruyama numerical scheme (Kloeden and Platen, 1992) at t n for a given θ and X 0 at t 0 with the interval [t n−1 , t n ] divided into M (large) subintervals of length g = (t n − t n−1 )/M . The multivariate normal density function, φ n,q (Y n,q |h ϕ (X (r) n ), θ) with expectation h ϕ (X (r) n ) and covariance Σ n (θ) given by equations (3.16 and 3.17) is where, n,q = Y n,q − h ϕ (X (r) n ) are the estimated errors at the time point t n for the q th experiment. Therefore, the simulated log-likelihood function based on equation (4.1) ( R (θ|Y ) = log L R (θ|Y )) is given by (4. 3) The parameter estimates are obtained by maximising the simulated loglikelihood function (4.3) with respect to θ.
Note that for S4, one does not need to simulate X n using Monte Carlo simulations but rather one solves the ODE, for example via a Matlab ode solver. To facilitate an efficient numerical search of the parameter space across a wide range of θ, the transformed parameter vector log 10 θ is estimated. Such estimates are robust since they are less prone to sampling error than when one directly works with untransformed parameters (Bland and Altman, 1996).

Model comparison
The likelihood ratio test (LRT) will be employed to compare two nested computational models in this paper (Madsen and Thyregod, 2010). Denoting the log-likelihood of the full and the reduced models by (θ Ha ) and (θ H 0 ), the likelihood ratio test statistic is given by Under certain regularity conditions (Wilks, 1962), LRT has an asymptotic (χ 2 ) chi-square distribution with (df ) degrees of freedom and df is the difference between the dimension of the vectors θ Ha and θ H 0 for large sample sizes.
Alternatively, for the non-nested models the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) given by AIC = 2(− (θ) + p), Anderson (2003) suggested the use of the corrected AIC (AIC C ) given by In fisheries studies, the AIC C has been shown to be accurate (Shono, 2000;Katsanevakis, 2006). Other related studies have also shown that the use of these information criteria is useful in comparing models (Shukor, 2014;Razzaghi et al., 2016;Mazerolle, 2006;Burnham et al., 2011). When comparing models, the one with the lower AIC, AIC C and BIC value is preferred.  As can be seen in Figure 3 (all panels), over the first 4 hours, there is a rapid decrease in the EE2-conjugate concentration, after which the rate of reduction slows considerably. The EE2-conjugate is converted to liquid EE2-parent (free), which increases rapidly from 30 mins to 8 hours and then slowly declines. The reason for zero concentrations being observed at t = 10 mins and t = 30 mins is that the measured concentrations may be below the analytical limit of detection, and so have been recorded as zero.

Experimental data
Whilst both of the above concentrations show transient behaviour, the EE2parent concentration in the solid phase (adsorbed or sequestered) does not, exhibiting only a steady slow growth over the 24 hours study period. Whilst the concentration in the solid phase was recorded as zero for times up to 2 hours, it may have been larger, but below the limit of detection.
The question this study aims to address is whether one can associate the steady slow growth trend for solid EE2-parent to adsorption or to absorption or to both. Since there is an adsorption -desorption process between liquid EE2-parent and sorbed EE2-parent, one may expect the profile of solid EE2parent to be similar to that of liquid EE2-parent but that is not the case here.

Model parameter estimation
The simulated log-likelihood estimates computed from equation (4.3) under hypotheses H1, H2 and H3 are presented in Table 2  When separate observational noise (3.17) is considered, the total number of parameters ranges from p = 8 for M2 under S4 to p = 13 for M0 under S1 for hypotheses H1 and H2 whilst hypothesis H3 have parameters ranges from p = 9 for M2 under S4 to p = 14 for M0 under S1. For this reason, the number of degrees of freedom, df , between any two nested models ranges from to 1 to 5. Thus, the critical values for at χ 2 df with df degrees of freedom at significance level of α = 0.05, χ 2 df (0.05) are 3. 8415, 5.9915, 7.8147, 9.4877 and 11.0704 respectively for df = 1, 2, 3, 4 and 5.
For each hypothesis, separate additive observational noise is tested against common additive noise, and it is found that a significance difference exists (i.e., separate additive observational noise is better than a common observational additive noise) across all models under hypotheses H1, H2 and H3.
The analysis of the models is further extended to examine the hypotheses H1, H2 and H3 with respect to separate additive observational noise.
Results from Table 2  -(4.6). Since C sor and C seq are usually indistinguishable, one might expect hypothesis H3 to hold, however results from the fit suggest that hypothesis H2 best describes the data. This outcome confirms previous findings on sequestration and in addition, quantified sequestered EE2 from the modelling perspective (Barret et al., 2011;Banihashemi and Droste, 2014).
The results of the simulated maximum likelihood estimation of the selected models of hypothesis H2 and S3 are shown in Table 3 and the dynamical simulated results with these parameters are depicted in Figure 4. The estimated value of the diffusion parameter or the system noise associated with the conjugate EE2 state variable (C con ), σ con is similar across all the three model structures considered (M0, M1 and M2). Statistically testing these mixed deterministic-stochastic models under S3 with purely deterministic models under S4, suggests that there are strong evidence for rejecting the purely deterministic models since the p-values corresponding to these tests are less than 10 −7 . Furthermore referring to Table 3, it can be observed that the difference between their estimated values of AIC, AIC C and BIC are small, thereby indicating that all these models adequately describe the experimental data. This is also evident in their simulated trajectories illustrated in Figure 4. Interestingly the model fits of M0, M1 and M2 under hypothesis H2 appear similar from the graphical point of view (see Figure 4). Table 3 Simulated maximum likelihood (SML) parameter estimate with 95% confidence intervals in parenthesis for hypothesis H2 under separate observational noise with respect to model structures (transformation pathways types) f (·) M0(k ab,f ree = 0, k ab,sor = 0), M1(k ab,f ree = 0, k ab,sor = 0) and M0(k ab,f ree = 0, k ab,sor = 0) System noise 3 (S3(σ f ree = 0, σ con = 0, σ sor = 0,  The EE2 SDE models M0, M1 and M2 under hypothesis H2 with respect to system noise type S3 capture the experimental data well over the entire 24 hours and, in addition, the empirical 95% confidence limits contain almost all the experimental data. From Figure 4, the uncertainty represented by the 95% empirical confidence limits especially for EE2-conjugate shows that the multiplicative noise in the SDE model accounts for influential factors that the deterministic ODE models cannot capture (not shown here).
This is evident in the parameter estimate of σ con given in Table 3. Therefore, the SDE model could account for certain inherent influential factors that are neglected by the ODE model (results under S4). However, it is worth noting that only one of the SDE models is subject to noise (dC con ) and the other models are deterministic in the chosen model. Figure S1 in the supporting information depicts histograms of the stochas- shows that there is no difference between M0, M1 and M2 as indicated by the likelihood ratio test and information criteria, and so M0, M1 and M2 are taken as plausible pathways. This suggests that more data are needed to distinguish these pathways. What can be said however is that, there is evidence of sequestration (Barret et al., 2011) and that the measured EE2-parent concentration in the solid phase, Y 3 corresponds to the EE2sequestered, C seq (which is hypothesis H2) for the experimental dataset (Gomes et al., 2009) used in this study. This is surprising, considering that in the lab and field, these two types of sorption to the solid phase (C sor and C seq ) are rarely defined separately (Joss et al., 2006). However, from the modelling point of view the solid (sorbed) parent compound (EE2), C sor equilibrates with the liquid parent compound C f ree so fast that the solid phase which is measured in the lab rather is the sequestered EE2 concentration in the solid phase (i.e., Y 3 = C seq ) as illustrated in Figure 4.
Thus, C seq is an important novel component in the model.

Perspective of ODE and SDE models
In general, the rich methodologies of both deterministic and stochastic modelling frameworks have been utilised to gain useful insights into the proposed models ( Figure 2) and to compare both approaches. The stochastic modelling results turn out to be significantly better (see Table 2) and further allows the distinction between the model system noise which is captured by the parameter σ con (accounting for inherent randomness of deconjugation of metabolites) and observational randomness accounted for in this study by separate additive noise.
Furthermore, the results of the parameter estimates obtained in this study show that not accounting for certain dynamics may lead to over-or under-estimation of model parameters which may contribute to variation in observed concentrations of emerging pollutants, as pointed out by authors (Liu et al., 2015;Barret et al., 2011). An example is the deconjugation kinetic rate coefficient, k dc obtained by Gomes et al. (2009) which is 0.32 h −1 and is less than the estimated value reported in Table 3 under hypothesis H2 where we find k dc = 0.42 h −1 . Differences in kinetic rate estimates in this paper to others are not surprising since the model accounts for sequestration dynamics, which to the best of the authors' knowledge has to date been neglected in the literature.

Significance
This present study investigates the mechanism of sorption of EE2 to the solid phase in an activated sludge process. In addition, the study examines the nature of the measured EE2 concentration in the solid phase as to whether it denotes the sorbed EE2 or sequestered EE2 or a combination of both from the modelling perspective. This study demonstrates and confirms sequestration, a phenomena known in literature yet rarely considered or distinguished beyond a pollutant associated with the solid phase (be that MLSS, sediment, sludge) (Barret et al., 2011;Yi et al., 2006;Yi and Harper, 2007). Results from the modelling strongly suggests that the experimentally measured EE2 concentration in the solid phase is the sequestered EE2, which is explained by hypothesis H2. In particular, for any BACs having high affinity for a solid phase like EE2 there may be similar behaviour because of similarities in their physiochemical properties.
This study shows the usefulness of mathematical modelling in wastewater treatment as the paper evaluates sequestration of EE2, and this will inform future efforts to develop and as well as apply models better able to fit and/or predict BAC behaviour and fate in the environment. The challenge, however, is that most data available from the literature does not capture the changes in the concentrations over time, where all of the three BAC states (e.g., liquid EE2-parent, solid EE2-parent and EE2-conjugate) are determined simultaneously. In the future, experiments of such nature may be encouraged to enable modelling as a route to better understand wastewater treatment process performance. Specifically, sequestered BAC concentrations will need to be measured as recommended in the literature (Barret et al., 2011).

Conclusions
Mathematical modelling of emerging pollutants in wastewater treatment offers a way of investigating phenomena and establishing which mechanisms affect the fate and behaviour of BACs in WWTPs. This study, has performed a proof-of-concept investigation of EE2 fate in an activated sludge pilot plant. In addition, the study examined the nature of the measured EE2 concentration in the solid phase via deterministic and stochastic modelling frameworks to gain insight into the fate and behaviour of emerging pollutants such as bioactive compounds, with a specific case study of EE2. Proposed models have been tested using the likelihood ratio test (LRT), Akaike Information Criterion (AIC), corrected AIC (AIC C ) and Bayesian Information Criterion (BIC) to examine the nature of the experimental measured EE2 concentration in the solid phase from a lab-scale batch experiment.
From the computational models developed in this paper, the proposed pathway of EE2 in wastewater treatment describes the experimental dataset well, thus confirming the existence of sequestration for the EE2 dataset. In addition, the results suggest that the experimental measured EE2 concentration in the solid phase is more likely to be the sequestered parent compound (C seq , hypothesis H2) rather than either the sorbed parent compound (C sor , hypothesis H1) or the sum of a fraction of the sorbed parent compound together with the sequestered parent compound (f sor C sor + C seq , hypothesis H3). Furthermore, the results based on information criteria in Section 4.2 show that purely deterministic models do not adequately fit the data and stochastic models of the process are required. However, not all components need to be treated stochastically, and this work indicates that the best fit was obtained when just one of the four variables (C sor ) was stochastic. This work shows that making the distinction between adsorbed and absorbed (sequestered) pollutants is beneficial. It can be hypothesized that fate studies involving pollutants that prefer the solid phase should evaluate sequestration phenomena in order to provide a more accurate description of the process.
The environmental consequence of this conclusion is that the uncertainty in EE2 or other BACs removal efficiency will be reduced, since wastewa-ter treatment can be further optimized with this insight and this will aid stakeholders in decision making.