Modeling the loading and fate of estrogens in wastewater treatment plants

Abstract Endocrine-disrupting compounds may produce infertility, nervous system disorders, and improper functioning of the immune system in humans and wildlife. Estrogens are classified as the most potent and common endocrine-disrupting compounds and the major point source for estrogen is municipal wastewater. Monitoring of estrogen is challenging, expensive, and intermittent; and therefore, the focus of this work is modeling estrone, 17β-estradiol, and 17α-ethynylestradiol concentrations from wastewater treatment plants in Calgary and Edmonton, Alberta, and Brandon, Manitoba. Demographic groups, excretion rates, population estimates, average daily flows, calculated estrogen transformation, calibration, calculated influent-to-effluent reduction percentages, and a treatment unit removal matrix are used to determine loading estimations of estrogen. Predicted average concentrations for EE2 and E2 in all the study sites exceed the threshold concentrations that could induce vitellogenin production by order of 13 and 2.3, respectively. The results demonstrate reasonable accuracy against previous measurements with r2 values ranging from 0.79 to 0.99 and RMSE values ranging from 0.5 to 9.4 ng/l and findings are consistent with concentrations reported in the literature. Upon further calibration with additional local data, the model may be used as a risk assessment analysis tool for these contaminants of concern.


PUBLIC INTEREST STATEMENT
Endocrine-disrupting compounds may produce infertility, nervous system disorders, and improper functioning of the immune system in humans and wildlife. Estrogens are classified as the most potent and common endocrine-disrupting compounds and the major point source for estrogen is municipal wastewater. In this paper, our focus was to model the fate and loading of the estrogens in the municipal waste water in Canadian context.

Introduction
Globally, there are growing health concerns regarding contaminants entering our water systems that provide us safe drinking water, irrigation water, food, and recreational opportunities (Wise, O'Brien, & Woodruff, 2011). Recently, occurrence of estrogens in aquatic environments have received considerable interest as they may lead to infertility, developmental disorders, disorders of the nervous system, and improper functioning of the immune system in humans and wildlife even at extremely low concentrations (National Institute of Environmental Health Sciences, 2015). Estrogens, estrone (E1), 17β estradiol (E2), and estriol (E3) are synthesized and excreted by humans (de Mes, Zeeman, & Lettinga, 2005). 17α estradiol (EE2) is a synthetic hormone which has molecular structure similar to that of E2, is the main estrogen used in oral contraceptives (OC), which is the most prescribed drug worldwide (Briciu, Kot-Wasik, & Namiesnik, 2009).
Concentrations as low as 0.5 ng/l of EE2 (Hansen et al., 1998;Purdom et al., 1994), 1 ng/l of E2 (Routledge, Parker, Odum, Ashby, & Sumpter, 1998), and 25 ng/L of E1 (Routledge et al., 1998) are reported to induce vitellogenin production in fish. In fact, recent studies have shown the presence of intersex fish downstream of wastewater treatment plants (WWTPs) at locations where estrogens are detected frequently and at higher concentrations among other endocrine disruptors (Evans, Jackson, Habibi, & Ikonomou, 2012). WWTPs are classified as a major point source for estrogens contaminants entering our waterways (Belfroid et al., 1999;Drewes, Heberer, Rauch, & Reddersen, 2003;Ternes, Kreckel, & Mueller, 1999). However, there are limited amounts of estrogen data collected since estrogens are currently unregulated pollutants (Umali, Pagsuyoin, & Parker, 2012). Furthermore, some cities in Canada are expanding and growing rapidly. The City of Calgary, for example, is expected to double to over 2.4 million residents by 2041 (Wright, 2014). The growing population will add increasing pressure on water quality and quantity in the region, including the Bow River (Robinson, Valeo, Ryan, Chu, & Iwanyshyn, 2009).
Measuring estrogen compounds in wastewater is challenging, expensive, and intermittent due to: the extremely low concentrations (pg-ng/L in water) of occurrence; complexity of the wastewater matrix; and the laborious analysis protocols (Johnson & Williams, 2004;Umali et al., 2012). Typical analytical procedure includes extraction (liquid-liquid or solid phase extraction), purification, hydrolysis, derivatization, and evaporation used in tandem with chromatography and mass spectrometry (Briciu et al., 2009;Umali et al., 2012). However, it has been previously acknowledged that monitoring estrogen concentrations in treated effluent is a requirement for accurate risk assessments (Jones, Voulvoulis, & Lester, 2001). Given that humans excrete estrogen on a daily basis (Johnson & Williams, 2004), and the major point source for estrogen are WWTPs (Belfroid et al., 1999;Drewes et al., 2003;Ternes et al., 1999), it should be practical to estimate the concentrations of estrogen that enter and leave WWTPs. Previous research by Johnson and Williams (2004) and others (Atkinson et al., 2012;Umali et al., 2012) have estimated estrogen concentrations using demographic group excretion rates, census statistics, flow rates, and WWTP operational parameters. Modeling estrogen is required in order to: (i) allow estimation of estrogen in avoiding high costs and time associated with laboratory analysis; (ii) better formulate a citywide risk estimate for these contaminants of concern; (iii) predict risks posed by individual WWTPs due to the presence of steroid estrogens in effluent in a city or region; and (iv) developing management strategies that focus on the contaminant(s) of most concern. As such, the objective of this research is to model the wastewater concentrations of E1, E2, and EE2 in order to predict the range of concentrations that can be expected to occur at WWTP's in Calgary and Edmonton, Alberta, and Brandon, Manitoba.

Study sites
In this study, estrogen release and its fate during wastewater treatment process was modeled in four selected WWTPs located within cities of Calgary and Edmonton, Alberta and in Brandon, Manitoba ( Figure 1). All the WWTPs chosen to simulate the model employ tertiary technologies; they have available data, and are in Western Canadian cities that are in relative close proximity to each other (Table 1a).   The Calgary wastewater flow data and historical population served by each WWTP was provided by the City of Calgary. The Brandon Manitoba WWTP parameters were obtained from Cicek, Londry, Oleszkiewicz, Wong, and Lee (2007), while the Edmonton WWTP parameters were obtained from Fernandez, Buchanan, and Ikonomou (2008) (Table 1b). The estrogen data used for comparison and validation to the simulated model were obtained from WWTPs located in Alberta and Manitoba ( Figure 1). As limited data were available for both Calgary facilities, supplemental measurements from similar size facilities were used for comparison and validation, i.e. Fish Creek WWTP measurements were supplemented with the Brandon WWTP measurements and the Bonnybrook WWTP measurements were supplemented with measurements from the Goldbar WWTP, Edmonton. The estrogen data collected to calculate the influent-to-effluent removal rates were obtained from various sources, mainly in Canada (Table 2); the estrogen removal matrix was calculated using a combination of estrogen data collected in Brandon, Manitoba (Cicek et al., 2007) and Wiesbaden, Germany (Andersen, Siegrist, Halling-sørensen, & Ternes, 2003) due to worldwide limited information. Data showing an influent-to-effluent concentration increase within WWTPs were not included in the model due to reported limitations in analytical methods preventing the detection of conjugated estrogen compounds commonly present in raw influent (Fernandez, Ikonomou, & Buchanan, 2007;Lee, Peart, & Svoboda, 2005).
The model calibration data-set included half of the EE2 influent and effluent measurements taken within the Brandon WWTP in 2003; the model validation data-set included the remaining EE2 measurements within the Brandon facility and also data from the Calgary facilities. Though estrogens are not detected in Calgary WWTPs during their wastewater analysis, it is assumed to be present but below the applicable detection limits. These concentrations below the detection limits were assumed to be zero and used to validate the calibrated model. The estrogen data collected from the WWTPs in Calgary and Brandon are compound-specific, while the estrogen data collected from the Edmonton WWTP in 2006 measures the total estrogenic activity or recombinant yeast assay (RYA) activity measured E2-eq data for the initial influent. Although, the RYA-measured E2-eq data correlate to the total estrogenic activity, E1 and E2 have been identified as the major contributors to RYA (Fernandez et al., 2008); however, since RYA measures the total estrogenic activity and lacks compound specificity, it was compared to the total estrogen concentration from the model or E1 + E2 + EE2. The detailed estrogen data collection methodology and analysis for the Brandon, Calgary, and Edmonton facilities are fully described in Cicek et al. (2007), Chen et al. (2006, and Fernandez et al. (2008), respectively.

Methods of data analysis
A schematic diagram illustrating the methods employed in this study is shown in Figure 3. The diagram consists of three main components: (1) revision of the demographic profile parameter ratios and calculation of the sewer biological transformation rate using local data, (2) influent and effluent model simulations and comparisons, and (3) calculated removal matrix and sensitivity analysis. Brief descriptions of the three major components are provided in the following paragraphs.

Demographic profile parameters
Demographic profiles for the model are divided into six categories due to differences in excretion rates of estrogen: menopausal females, menstruating females, females taking hormone replacement therapy (HRT), pregnant females, females on birth control pills (BCP), and males. The excretion rates for the six demographic profiles used in this study were taken from Johnson and Williams (2004) from two standard errors to give 95% confidence limits; however, the group sizes have been revised based on local data (Table 3). In addition, the male and female parameter definitions were revised to exclude prepubescent individuals, i.e. 13 years of age and under, due to extremely low estrogen excretion concentrations (Jorgensen, Keiding, & Skakkebaek, 1991). The parameter revisions should improve model accuracy over previous simulations completed by Johnson and Williams (2004), where excretion rates were assumed to be constant over the entire female population and did not exclude prepubescent individuals.

Biological transformation rate of estrogens during transit in sewers
Researchers have reported that E2 can readily biodegrade to E1 during sewer transit due to the presence of native biofilms in sewers that have metabolic capabilities associated with biodegradation (D'ascenzo et al., 2003;de Mes, 2007;Johnson & Williams, 2004;Ternes et al., 1999). In contrast, EE2 degradation occurs at much slower rates and in some cases no degradation occurs due to the presence of the ethinyl group that's responsible for delaying enzyme expression, substrate-receptor binding, and metabolism (de Mes, 2007;Racz & Goel, 2009). In fact, previous observations in activated sludge microcosms have shown no transformation of E1 & EE2 during sewer transit (Johnson & Williams, 2004). Therefore, only the transformation of E2 to E1 will be considered for this study. The estimated sewer biological transformation rate of E2 to E1 (i.e. R T ) for this study was calculated by (1) running the model for the Brandon WWTP due to available compound-specific influent data, (2) assessing whether the model agreed with the data when analyzing E1 + E2 together against E1 + E2 measurements, and (3) if in agreement, calculate the difference between the modeled E2 and the measured E2. The difference between the modeled E2 and the measured E2 would therefore constitute the transformation rate of E2 to E1. The preliminary model simulation results show agreement with the measurements as all of the daily mean measurements for E1 + E2 fell within the predicted intervals, and the estimated biological sewer transformation rate of E2 to E1 was therefore calculated to be 27% based on the percent difference equation: where E 2,a and E 2,b represents the mean modeled value and mean measurements, respectively. (1)

EE2 Model Calibration
The Johnson and Williams (2004) model will form the basis used to determine the total amount of estrogen arriving at the WWTP: where T E (μg/l) is the total estrogen influent loading concentration, T E1 (μg/l) represents the total E1 influent loading concentration, T E2 (μg/l) represents the total E2 influent loading concentration, and T EE2 (μg/l) represents the total EE2 loading concentration.
where R T represents the percentage of E2 removed in the sewer system, P G,i is the number of people within the demographic group of the WWTP served population, U Ei (μg/d) is the estrogen excretion rate in urine per individual of a particular demographic profile, F Ei (μg/d) is the estrogen excretion rate in feces per individual of a particular demographic profile, and E c is the concentration of E1 transformed from E2 during the sewer transit (i.e. the equivalent concentration of E2 removed in the sewer system represented as R T ), and Q is the WWTP average daily flow rate (l/d). R T is estimated to be 27% based on the results of the Brandon model preliminary simulation. E c for Equation (4) will be the amount of degraded E2 (i.e. 27%). The excretion value rates used account for the free, For E1:

Model effluent simulation
The following equation describes the removal efficiency calculation based on the raw measured influent and effluent data as referenced in Table 2: where R E is the percentage of the WWTP influent load removed, E infl is the influent load (μg/l), and E Effl is the effluent load (μg/l).
The mean percent estrogen influent-to-effluent removal values were calculated with 95% confidence limits using the R E results calculated from Equation (7) above. For example, to determine the highest predicted effluent concentration, with 95% confidence limits, the upper standard error value for EE2 excretion concentration and the lowest predicted removal performance are used; to determine the lowest predicted effluent concentration, the lower standard error value for EE2 excretion concentration and the highest predicted removal performance are used ( Table 2). The following equation describes how the 95% confidence limits were calculated using the mean percent estrogen influent-to-effluent removal values: where CL 95 is the 95% confidence limit, X is the mean percent estrogen influent-to-effluent removal values, SD is the standard deviation, and n is the sample size.
The following equation describes the final effluent calculation based on the estimated influent concentration and calculated removal percentage efficiency: where F EEffl (μg/l) is the final estimated effluent concentration value for a particular estrogen compound.

Calculated removal matrix and sensitivity analysis
The WWTP concentrations estimated to be found after each individual treatment process was calculated using data from the Brandon and Wiesbaden WWTPs. Since some components of the Brandon and Wiesbaden treatment process differs from Bonny Brook and Fish Creek, only the data obtained from relevant process units was included in the calculations. For example, the Brandon WWTP employs a high-intensity UV disinfection unit similar to Fish Creek and Bonnybrook, whereas the Wiesbaden WWTP does not employ UV disinfection. As a result, UV disinfection data from the Brandon facility were used. In contrast, the Brandon WWTP employs a gravity vortex grit removal unit as part of its primary treatment, whereas Wiesbaden employs a primary clarifier (i.e. open-air tanks used for settling and skimming) similar to Bonnybrook and Fish Creek. As a result, data from the Wiesbaden WWTP Primary treatment unit were used verses the vortex grit removal process found in Brandon. The percent difference was calculated using Equation (1) as noted above.
The sensitivity index (SI), developed by Hoffman and Gardner (1983), accounted for the possible values when determining parameter sensitivity: where max(Pi) and min (Pi) are maximum and minimum output values, respectively, from the range of input values used. The SI for each demographic profile and associated excretion rate variables were completed within the upper and lower range of the given value (see in Table 4).
The model was validated based on the following criteria: (i) the 95% confidence intervals for the estrogen influent and effluent predicted by the model for a given WWTP should cover 90% of all daily mean estrogen measurements-a similar approach used by Ram and Gillett (1993); (ii) WWTP's with limited estrogen data was validated using supplemental data from similar WWTP's where 90% of all daily mean estrogen measurements lie within the 95% confidence intervals; and (iii) a regression analysis between measurements and predictions. Figure 4 shows a comparison between the modeled WWTP influent and effluent verses measurements. The model produced good results as all of the E1 and E2 influent and effluent mean measurements per day are within the upper and lower confidence limits of the model. Discrepancies between the predicted and measured concentrations for EE2 may be due to underestimating the number of people taking OC and deviations from the average daily flow (Johnson, Belfroid, & Di Corcia, 2000). Only through extensive calibration, does the model allow for all of the mean measurements of EE2 to fall within the upper and lower limits. Model calibration was achieved by adjusting the upper and lower limits within the Brandon WWTP model to a point that included all of the measurements for half of the Brandon data-set while running the model against the remaining samples. In this case, the EE2 influent upper limits increased by factor of 6.4 and the lower limits decreased by a factor of 1.4; for EE2 effluent, the upper and lower removal percentage decreased to 62% and 0%, respectively. The calibration model for EE2 was also partially validated using effluent data from  the Calgary facilities. For the Goldbar facility, 79% of the mean daily influent estrogen simulated values lie within the 95% confidence intervals for the model. The individual influent measurements that fell outside of the confidence limits, could be explained by the presence of other potential compounds producing estrogenic effects as well as competing estrogenic and anti-estrogenic substances, both of which E2-eq estrogen measurements are known to account for. Overall, the model shows agreement with the measurements, and the model results are also consistent with concentrations reported in the literature (Baronti et al., 2000;Fernandez et al., 2007;Lee et al., 2005;Lishman et al., 2006;Servos et al., 2005).

Results and discussion
As mentioned previously, concentrations to induce vitellogenin production in fish have been reported as low as 0.5 ng/l of EE2 (Hansen et al., 1998;Purdom et al., 1994), 1 ng/l of E2 (Routledge et al., 1998), and 25 ng/L of E1 (Routledge et al., 1998). The model concentrations for EE2 and E2 averaged in all the study sites exceed the threshold concentrations by order of 13 and 2.3, respectively. These predicted concentrations suggest potential risk to fish may exist as vitellogenin production in fish could be induced. Furthermore, cumulative impacts of WWTPs, agriculture, landfills, industry, and humans also may contribute to the potential risk. However, in contrast, seasonal variations in temperature, hydrology, contaminant dispersion, and river dilution ratios could dramatically increase or decrease concentrations and need to be taken into consideration to accurately assess risk. Figure 5 shows the relationship between measured E2-eq. data and total estrogen modeled values for the Edmonton Goldbar facility; Figure 6 shows the relationship between compound specific data for E1, E2, and EE2. The model produced good results with influent and effluent r 2 values ranging from 0.79 to 0.99 and RMSE values ranging from 0.5 to 9.4. However, modeling exercises in the past has been confounded by conflicting results. In one modeling study using data from five Italian and three Dutch WWTPs, Johnson et al. (2000) produced reasonable influent predictions for E1 (r 2 value of 0.50) and E2 (r 2 value of 0.47), while EE2 (r 2 value of 0.149) proved to be less accurate. In another study, Johnson and Williams (2004) modeled estrogen for six Roman and one German WWTP, and produced reasonably good results for predicting influent concentrations with r 2 values of 0.70 (p < 0.001) for E1, 0.66 (p < 0.001) for E2, and 0.53 (p = 0.06) for EE2. In addition, effluent predictions for E2 and EE2 all fell within the ranges of observed values, while E1 proved difficult to predict. Umali et al. (2012) modeled E2 influent loading for WWTPs in Bethlehem and Allentown, and produced good results that varied 1.2-2.5% from observed values for Bethlehem and Allentown, respectively. In contrast, Atkinson et al. (2012) found that measured E1 and E2 concentrations in wastewater influent for WWTPs in Ottawa and Cornwall were higher than predicted estimates-by a factor of 2 for E1.

(a) Influent (b) Effluent
As discussed previously, deviations between the model and measurements for the Goldbar facility could be explained by the presence of other potential compounds producing estrogenic effects as well as competing estrogenic and anti-estrogenic substances, both of which E2-eq estrogen measurements are known to account for. Also, reported limitations in analytical methods are known to prevent the detection of conjugated estrogen compounds commonly present in raw influent (Fernandez et al., 2007;Lee et al., 2005), whereas the model accounts for unconjugated estrogen compounds. Deviations between the model and measurements for the Brandon facility could be explained by sampling times during the day which may capture peak estrogen loading, deviations in flow data, and grab samples that do not represent daily averages. Note that limited available data for the Calgary facilities prevented the completion of a regression analysis, and suggested that additional data would be required to better corroborate the Calgary models.

Calculated removal matrix and sensitivity analysis
The results of the removal matrix simulation show that the amount of estrogen is subsequently reduced within each successive treatment unit within the WWTP, with the activated sludge bioreactor showing the greatest reduction (Table 5 and Figure 7). The results are consistent with the literature. For example, Matsui et al. (2000) reported that the estrogenic activity measured by yeast estrogen screening decreased after each processing unit within the WWTP, and denitrification within the activated sludge treatment process showed the greatest decrease. Similarly, in Germany, Andersen et al. (2003) reported lower concentrations of estrogen with each successive treatment unit within the Wiesbaden WWTP, and the greatest reduction occurred at the denitrification stage. It is explained that denitrification and dilution-with the return sludge from the secondary clarifier and the internal recirculation containing little estrogen from the last nitrification tank-contributed to the reduced concentrations of estrogen (Andersen et al., 2003). The reductions also suggest that de-nitrification and aerobic biological degradation appear to play a role in reducing the amount of E1/E2 and EE2, respectively. The reduced concentrations may also suggest slow sorption kinetics and no equilibrium between the absorbed and dissolved estrogens. The final effluent concentrations of the removal matrix are also supported by the model as all of the mean concentrations fall within the upper and lower limits of the model for the same year.  Finally, a sensitivity analysis was completed by calculating the estrogen influent loading concentration for E1, E2, and EE2 while varying each individual variable to its upper and lower range (Table 4); these calculations represent the maximum and minimum output values max (Pi) and min (Pi) used in the SI equation (Equation (10)) to determine the output % difference. The results for E1 can be expressed as follows: pregnant excretion rate > pregnant females > menstruating female excretion rate > daily treated flow of wastewater > females on HRT > menstruating female excretion rate > male excretion rate > HRT excretion rates > menstruating females > menopausal females > males; the results for E2 can be expressed as follows: pregnant females > pregnant excretion rate > daily treated flow of wastewater > females on HRT > menstruating female excretion rate > male excretion rate > menopausal female excretion rate > HRT excretion rates > menstruating females > menopausal females > males; and the results for EE2 can be expressed as females taking BCP > daily treated flow of wastewater > EE2 excretion rate (Figure 8). The results of E1 and E2 are generally consistent and the noted deviations can be explained by variances in the excretion limits. Overall, the percentage of females taking the BCP shows the greatest sensitivity variance of approximately 40%; and therefore, is a variable of importance when modeling EE2.

Conclusions
In this paper, a simple model was employed to determine whether an accurate range of estrogen can be predicted to arrive and leave a WWTP in Western Canada. This study employed (i) estrogen excretion rates and revised parameter ratios (ii), model simulations using a calculated sewer biological transformation rate and influent-to-effluent removal rates, and (iii) a calculated removal matrix and sensitivity analysis. The analysis revealed reasonably good results for predicting E1 and E2 values and calibration was required to produce good results for EE2. Rate of excretion of E1 during pregnancy, population of pregnant females and population taking OC are found to be the most sensitive parameters that influence the model predictions. Improvement in model predicted E1, E2, and EE2 concentrations suggest that the mathematical model presented here could be used to assist with formulating a citywide risk estimate for these contaminants of concern.