A Closer Look on Spatiotemporal Variations of Dissolved Oxygen in Waste Stabilization Ponds Using Mixed Models

Abstract

Dissolved oxygen is an essential controlling factor in the performance of facultative and maturation ponds since both take many advantages of algal photosynthetic oxygenation. The rate of this photosynthesis strongly depends on the time during the day and the location in a pond system, whose roles have been overlooked in previous guidelines of pond operation and maintenance (O&M). To elucidate these influences, a linear mixed effect model (LMM) was built on the data collected from three intensive sampling campaigns in a waste stabilization pond in Cuenca, Ecuador. Within two parallel lines of facultative and maturation ponds, nine locations were sampled at two depths in each pond. In general, the output of the mixed model indicated high spatial autocorrelations of data and wide spatiotemporal variations of the oxygen level among and within the ponds. Particularly, different ponds showed different patterns of oxygen dynamics, which were associated with many factors including flow behavior, sludge accumulation, algal distribution, influent fluctuation, and pond function. Moreover, a substantial temporal change in the oxygen level between day and night, from zero to above 20 mg O₂·L⁻¹, was observed. Algal photosynthetic activity appeared to be the main reason for these variations in the model, as it was facilitated by intensive solar radiation at high altitude. Since these diurnal and spatial patterns can supply a large amount of useful information on pond performance, insightful recommendations on dissolved oxygen (DO) monitoring and regulations were delivered. More importantly, as a mixed model showed high predictive performance, i.e., high goodness-of-fit (R² of 0.94), low values of mean absolute error, we recommended this advanced statistical technique as an effective tool for dealing with high autocorrelation of data in pond systems.

1. Introduction

Waste stabilization ponds (WSPs) have increasingly received attention since these shallow lagoons offer a natural biological purification of wastewater with low cost and minimal operation and maintenance requirements [1]. In fact, thousands of its applications currently serve millions of people in many countries across the globe. For example, Dandora WSP in Kenya, the biggest pond treatment system in Africa, serves approximately one million inhabitants or the wastewater from a population of 1.6 million is treated by Western WSP at Werribee in Melbourne, Australia [2]. A distinctive factor of WSPs which allows differentiating them from conventional wastewater treatment plants (WWTPs), is the involvement of algal photosynthesis. During the day, the algal photosynthetic process generates oxygen for aerobic heterotrophs to mineralize organic matters, which, in turn, produce CO₂ for the growth of algae [3]. Taking advantage of this natural oxygenation, pond treatment systems reduce operational costs and constrain potential risks from the emission of volatile organic compounds by avoiding mechanical aerations [4]. On the other hand, during the night or under light-limited conditions, such as in cloudy days or at certain water depths, instead of oxygenation, algae respire and thus consume oxygen and release CO₂ [3]. In short, the metabolism of algae, which is strongly dependent on spatiotemporal properties and meteorological conditions, can cause a wide variation of the oxygen level in WSPs [5]. Therefore, it is not an easy task for pond engineers to have a proper regulation of oxygen level in which respiratory oxygen required by aerobic bacteria is met by algal photosynthetic oxygen without any additional mechanical aerations.

Although pond technology has been developed over decades, the number of models serving for a better understanding of oxygen dynamics in pond systems remains small. To the authors’ knowledge, there are only two studies, i.e., Kayombo et al. [6] and Banks et al. [7], which applied mathematical models to investigate the oxygen balance in facultative ponds (FPs). Including only algal photosynthesis, the model of Kayombo et al. [6] considered four driving forces of oxygen variation in pond systems, i.e., light intensity, pH, temperature, and CO₂. This model suggested that 99% of oxygen production was from the algal photosynthesis while the inflow from primary FP brought 1% left. Banks et al. [7] advanced their model by adding aerobic bacterial assimilation of organic matter whose rate was strongly affected by temperature. Although both studies considered the effect of climatic factors, i.e., light intensity and air temperature, on oxygen balance in pond systems, the interactions between these climatic factors and temporal characteristics were not taken into account. Particularly, even though the hourly variations of light intensity and water temperature were clearly depicted in both studies, only the daily average values were applied in the models instead. In addition, wind mixing as the second mechanism of oxygenation was also neglected in these models. This mass diffusion from the atmosphere was considered as a predominant influencing factor on WSP performance in Li et al. [8].

Oxygen dynamics in pond systems is in terms of not only time but also space. The spatial variation of oxygen is associated with the change in algal community composition and distribution between different ponds and different locations within each individual pond [9]. Pham et al. [10] observed higher values of algal abundance, richness and diversity in maturation ponds (MPs) compared to FPs and lower biovolumes of motile algal species in the inlet compared to the outlet of Ucubamba WSP system in Cuenca (Ecuador). Furthermore, as a result of light attenuation, the algal photosynthesis in FPs can locate only at 20–30 cm from the water surface while that value is around 60 cm for clear and less turbid MPs [11]. In short, these different distributions of algae and its stratification in different depths create a dynamic spatial pattern of dissolved oxygen (DO) in pond systems.

Therefore, our main objective is to investigate the spatial and temporal effects on oxygen dynamic in the WSPs. To this end, the first application of linear mixed effect models (LMMs) in WSPs was implemented. Thanks to its ability to analyze clustered longitudinal data and repeated measures, the spatial and temporal autocorrelations of data can be taken into account [12]. This model was fitted on the data collected from three meticulous sampling campaigns in Ucubamba WSP in Cuenca (Ecuador). Especially noteworthy is that this pond system was located at high altitude, i.e., 2400 m a.s.l., where climatic conditions are relatively severe, i.e., strong solar radiation, low oxygen pressure, and low air temperature with high variation [13]. Therefore, it is expected that the combination effect between the spatiotemporal variations and these extreme meteorological conditions can generate a very dynamic pattern of oxygen level in this WSP. Subsequently, the predictive performance of the mixed model was evaluated via leave-one-out cross validation (LOOCV), through which the detailed description of oxygen dynamics at different depths and daytimes was illustrated. More importantly, considering these findings, together with the fact that DO control in WSPs is overlooked in the previous guidelines of pond operation and maintenance (O&M), insightful recommendations for DO monitoring and control in WSPs were drawn.

2. Materials and Methods

2.1. Study Area

Ucubamba WSP, the biggest wastewater treatment plant in Ecuador, is operated to purify the domestic wastewater of more than 400,000 inhabitants in Cuenca. The pond system is located at the altitude of 2400 m a.s.l. where annual average temperature is 14 °C and average rainfall is about 780 mm per year. Subsequent to pre-treatment steps including screening and grit chamber, the pond system is divided into two identical flow lines, each of which contains an aerated pond (AP), a facultative pond (FP) and a maturation pond (MP) (Figure 1). As a primary treatment, the APs are used as an alternative for anaerobic ponds to remove organic matter and nutrients [14]. These two ponds have a depth of 4.5 m with aerators to supply enough oxygen for oxygenation but not to keep the biomass and influent materials in suspension, hence, organic matters can still be decomposed anaerobically at the bottom of the ponds [15]. With a depth of 1.7 m, FPs decompose the organic matters mainly under aeration conditions as a result of the oxygenation from the microalgal photosynthetic activity. MPs (1.5 m depth) further polish the wastewater, especially in terms of pathogen removal.

2.2. Sampling Campaigns

Three sampling campaigns were conducted on 25/26 July 2013, 14/15 August 2013 and 26/27 August 2013. Each sampling campaign lasted for two days as the sampling of one line required the period of one day, i.e., from 8:00 a.m. until 6:00 p.m. This course of time covers the whole period of daylight in Cuenca, ensuring the investigation of temporal effects on algal metabolism, hence, on the oxygen variation. Each pond was divided into six sections along the longitudinal direction and four sections breadthways, which created three zones: influent (location 1, 2 and 3), middle (location 4, 5 and 6), and effluent (location 7, 8 and 9) (Figure 1). At two different depths, 30 cm below the water surface and 15 cm above the sediment layer, temperature (°C), pH (−), chlorophyll a (µg Chl a·L⁻¹) and DO (mg O₂·L⁻¹) were determined by two manual multi-probes, YSI 6600 V2 and YSI 6920 V1. These probes were carefully calibrated every three days by following their manual in order to ensure their accuracy. At the same time, mixed samples of each zone were analyzed at two different depths for biochemical oxygen demand (BOD₅, mg O₂·L⁻¹), chemical oxygen demand (COD, mg O₂·L⁻¹), total Kjeldahl nitrogen (TKN, mg N·L⁻¹), total phosphorus (TP, mg P·L⁻¹), and total solids (TS, mg·L⁻¹) using American Public Health Association methods [17]. Due to the sludge accumulation, the samples at the bottom of location 1 and 2 of the FPs could not be collected. Meteorological data, including air temperature (°C), solar radiation (W·m⁻²) and wind speed (m·s⁻¹), were obtained from the meteorological station of CELEC Hidropaute, located 600 m away from the WSP.

2.3. Kruskal–Wallis and Bonferroni Correction

Before applying the mixed model, Kruskal–Wallis tests followed by Bonferroni–Dunn test were applied for multiple comparisons of oxygen between different sampling campaigns, different locations within a pond and among the ponds. Unlike parametric tests, such as one-way ANOVA test, the Kruskal–Wallis test, a non-parametric statistical tool, does not require an assumption of normal distribution of residual. To avoid a type I error, Bonferroni correction is widely applied for correcting the p-values in multiple comparison tests [18]. These tests were carried out using “dunn.test” package [19] in R software Version 3.0.2 [20]. The p-value was considered significant at 0.05/n with n as the number of hypotheses being tested in multiple comparisons.

2.4. Model Selection

One of the main objectives of our research is to investigate the effects of spatiotemporal characteristics and their interactions with meteorological conditions on the variation of oxygen level within the WSP. To this end, LMM, as an advanced technique for statistical modeling, was executed in R [20] using the lme function in the nlme package [21]. Not only taking into account fixed effect as linear regression models, LMM are comprised of both fixed effects and random effects, which can take into account the spatiotemporal autocorrelations of data [22]. The determination of fixed-effect variables was based on the mass balance of oxygen within the ponds. While the main oxygen sources in the WSP system were photosynthesis and the direct exchange of atmospheric oxygen through the air/water interface, oxygen consumption was mostly done by aerobic bacteria for mineralizing organic matter and nitrification process [23]. Particularly in the model, chlorophyll a (µg·L⁻¹) and solar radiation (W·m⁻²) characterized the photosynthetic activity, while wind speed (m·s⁻¹) and air temperature (°C) represented the oxygen exchange processes. BOD₅ and COD represented the bacterial mineralization whilst TKN and TP were nutrients for the growth of algae and bacteria and nitrification process. The spatial and temporal variation of these variables in ponds were also reported in previous studies, e.g., McLaughlin et al. [24] and Guo et al. [25]. Moreover, we modelled the effects of depth and daytime as a logarithm function and a quadratic function, respectively, based on their observed patterns in the studies of Kayombo et al. [6] and Tadesse et al. [26]. Most importantly, the interactions between daytime and the three meteorological parameters were also simulated in the LMM.

Regarding random effect, pond and sampling-campaign parameters were included to account for the spatial and temporal autocorrelation between samples, creating a three-level hierarchical mixed model. More specifically, the unit of analysis, DO concentration (level 1), is nested within pond (level 2), which is in turn nested within sampling campaign (level 3). The detail of this three-level mixed model is demonstrated in Equation (1) and Figure 2:

$$\begin{array}{l}{\mathrm{LogDO}}_{\mathrm{ijk}}{ = \mathsf{\beta }}_0{ + \mathsf{\beta }}_1{\times \mathrm{BOD}}_{\mathrm{ijk}}{+ \mathsf{\beta }}_2{ \times \mathrm{COD}}_{\mathrm{ijk}}{ + \mathsf{\beta }}_3{ \times \mathrm{TS}}_{\mathrm{ijk}}{ + \mathsf{\beta }}_4{ \times \mathrm{TN}}_{\mathrm{ijk}}{+ \mathsf{\beta }}_5{ \times \mathrm{TP}}_{\mathrm{ijk}}{ + \mathsf{\beta }}_6{ \times \mathrm{pH}}_{\mathrm{ijk}}{ + \mathsf{\beta }}_7{ \times \mathrm{Chl}}_{\mathrm{ijk}} +\\ {\mathsf{\beta }}_8 \times {\text{Solar.rad}}_{\mathrm{ijk}}{ + \mathsf{\beta }}_9 \times {\text{Wind.speed}}_{\mathrm{ijk}}{ + \mathsf{\beta }}_{10} \times {\text{Air.temp}}_{\mathrm{ijk}}{ + \mathsf{\beta }}_{11}{ \times \mathrm{LogDepth}}_{\mathrm{ijk}} +{\mathsf{\beta }}_{12}{ \times \mathrm{Daytime}}_{\mathrm{ijk}}\\ + {\mathsf{\beta }}_{13}{ \times \mathrm{Daytime}}_{\mathrm{ijk}}^2{ + \mathrm{Daytime}}_{\mathrm{ijk}}{ \times (\mathsf{\beta }}_{14} \times {\text{Solar.rad}}_{\mathrm{ijk}}{ + \mathsf{\beta }}_{16} \times {\text{Wind.speed}}_{\mathrm{ijk}} +{\mathsf{\beta }}_{18} \times {\text{Air.temp}}_{\mathrm{ijk}})\\ {+\mathrm{Daytime}}_{\mathrm{ijk}}^2{ \times (\mathsf{\beta }}_{15} \times {\text{Solar.rad}}_{\mathrm{ijk}}{ + \mathsf{\beta }}_{17} \times {\text{Wind.speed}}_{\mathrm{ijk}}{ + \mathsf{\beta }}_{19} \times {\text{Air.temp}}_{\mathrm{ijk}}) {+ \mathrm{a}}_{\mathrm{k}}{ + \mathrm{a}}_{\mathrm{jk}}{ + \mathsf{\epsilon }}_{\mathrm{ijk}}\end{array}$$

(1)

where:

${\mathrm{DO}}_{\mathrm{ijk}}$: The concentration of DO of observation i within pond j collected at sampling campaign k. i: 1–16 for the FPs, 1–18 for the MPs, j: 1–4; k: 1–3,

${\mathsf{\beta }}_1- {\mathsf{\beta }}_{13}$: Fixed effects of the 12 variables,

${\mathsf{\beta }}_{14}- {\mathsf{\beta }}_{19}$: Interaction between daytime and the meteorological variables,

${\mathrm{a}}_{\mathrm{k}}$: Random effect associated with the intercept for sampling campaign k. ${\mathrm{a}}_{\mathrm{k}}$ ~ N(0,${\sigma }_{sc}^2$),

${\mathrm{a}}_{\mathrm{jk}}$: Random effect associated with the intercept for pond j within sampling campaign k. ${\mathrm{a}}_{\mathrm{jk}}$ ~ N(0,${\sigma }_p^2$),

${\mathsf{\epsilon }}_{\mathrm{ijk}}$: Residual. ${\mathsf{\epsilon }}_{\mathrm{ijk}}$ ~ N(0,${\sigma }^2$).

To assess the assumptions of the mixed model, such as outliers, multicollinearity, and homogeneity, data exploration was performed following the guidelines of Zuur et al. [27]. Prior to statistical analyses, we assessed the followed assumptions: (1) outliers via the means of Cleveland dotplots (Figure S1), (2) multicollinearity using pairwise scatter-plots comparing the correlation coefficients among covariates (Figure S2), and (3) homogeneity via the residuals of the fitted model (Figures S3 and S4). After removing outliers, the remaining DO concentrations were log10 transformed to stabilize the variance for statistical analyses and then the model was fitted with the dataset. We evaluated the goodness-of-fit of the model via conditional R² for both fixed and random effects [28]. For residual diagnostics, normality and homogeneity were tested via the QQ plot and residuals vs. fitted values plot (see Figures S5–S7).

2.5. Model Evaluation

The predictive performance of mixed model was evaluated by the mean absolute error (MAE) using leave-one-out cross validation (LOOCV). In particular, the MAE, representing the mean deviation between observed values and predicted ones, was calculated as follows (Equation (2)) [29]:

$$\mathrm{MAE}=\frac{1}{\mathrm{n}}{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{O}}_{\mathrm{i}}-{\mathrm{P}}_{\mathrm{i}}\right|,$$

(2)

where O_i is the observed DO in sample i, and P_i is the corresponding prediction based on the mixed model fitted with the full dataset of n samples but without sample i (LOOCV). MAE was chosen over root mean square errors (RMSE) since MAE was concluded as the most natural measure of average error in contrast to the inconsistent functional relationship between RMSE and average error, which might lead to confused interpretations [30,31]. Moreover, the bias and consistency of model prediction were evaluated by regressing observed vs. predicted oxygen concentrations [32].

2.6. Intraclass Correlation Coefficient (ICC)

Autocorrelations in time and space appear when the values of data sampled at the same time and location exhibit more similar patterns than those at different sampling times or further apart. Without considering spatial and temporal autocorrelations, the linear regression model can violate the assumption of independently and identically distributed random variables and draw incorrect conclusions [22]. On the other hand, mixed models with random effects can represent the impact of these autocorrelations by the mean of intraclass correlation coefficient (ICC), which is a measure describing the homogeneity of the observed oxygen concentrations within given clusters, i.e., pond and sampling campaign [33]. ICC is determined as a function of the variance components in a mixed model. For example, sampling-campaign-level intraclass correlation coefficient, ICC_sc, was calculated by dividing the variance of the random sampling-campaign effects (${\sigma }_{sc}^2$) by the total random variation. The latter consisted of ${\sigma }_{sc}^2$, the variance of the random effects associated with ponds nested within sampling campaign (${\sigma }_p^2$) and the variance of residual (${\sigma }^2)$ (Equation (3)):

$${\mathrm{ICC}}_{\mathrm{sc}}=\frac{{\sigma }_{sc}^2}{{\sigma }_{sc}^2+{\sigma }_p^2+{\sigma }^2}.$$

(3)

The value of ICC_sc is high when the total random variation is dominated by ${\sigma }_{sc}^2$, meaning that the oxygen concentrations measured among different sampling campaigns tended to widely vary while these values among different ponds within a sampling campaign are relatively homogenous. The pond correlation coefficient, ICC_pond, was calculated as the proportion of the variance of the random effects, ${\sigma }_{sc}^2+{\sigma }_p^2$, to the total random variation (Equation (4)):

$${\mathrm{ICC}}_{\mathrm{pond}}=\frac{{\sigma }_{sc}^2+{\sigma }_p^2}{{\sigma }_{sc}^2+{\sigma }_p^2+{\sigma }^2}.$$

(4)

The pond-level ICC is high if there is little variation in the oxygen level within the same pond relative to the total random variation (${\sigma }^2$ is low).

3. Result

3.1. Spatial Variation of Dissolved Oxygen

The spatial variations of oxygen level at different ponds and depths are demonstrated in Figure 3. A wide variation of DO was found among four ponds, which corresponded to low p-values of Bonferroni–Dunn tests (p-values < 0.005), except for the comparison between two ponds FP2 and MP2 (p-value = 0.8526). Indeed, there was a different behavior between the two flow lines. Particularly, FP1 contained higher concentrations of DO than its counterpart at the top line, but these values of its consecutive pond (MP1) significantly declined and were lowered in MP2. In fact, from the outlet part of FP1 to MP1 inlet, DO values near the water surface dropped about 70%, i.e., from above 10 mg O₂·L⁻¹ to around 3 mg O₂·L⁻¹, while the oxygen level remained similar between two ponds in the upper line. This trend also occurred at the bottom layers of the ponds. Between two depths, higher values of both concentration and variation were found close to the surface as being supported by a low p-value of Kruskal–Wallis test (p-value < 0.0001). Much less expected is an extremely high value of oxygen concentration at the bottom of FP1, around 17 mg O₂·L⁻¹, which was observed during the last two hours of the afternoon in the first sampling campaign. These extreme values caused very high deviation of the oxygen concentration at the bottom of FP1.

For a further investigation, DO concentrations and variations at the nine locations of each pond in the system are illustrated in Figure 4. Via this bubble plot, the variations of oxygen level within a pond and between two depths are evidently showed. As such, heterogeneous oxygen concentrations were found at different zones across the water surface of the FPs. For example, at the surface of FP1, the highest concentrations were located in the middle area with around 4 mg O₂·L⁻¹ higher than those values in the influent and effluent area. Likewise, we also observed higher concentrations and fluctuations at FP2. On the other hand, the oxygen level in the MPs remained homogeneously, around 5 mg O₂·L⁻¹ at MP2 surface and 3 mg O₂·L⁻¹ at MP1.

3.2. Model Selection

Prior to the statistical modeling, we used pairwise scatter-plots to compare the correlation coefficients among covariates to the threshold of 0.7 as a suitable indicator for the severe distortion effects on model estimation caused by collinearity [34]. The statistical analysis showed that the correlation coefficients among three parameters, i.e., daytime, wind speed, and air temperature, were larger than the threshold (see Figure S2). Hence, we dropped the two meteorological parameters as daytime parameters can be measured with the least effort and cost [27]. Likewise, BOD, COD, and TS also showed high multicolinearity; thus, we removed BOD and TS from the model with the same reason.

In the next step of data exploration, the assumption on homoscedasticity of residual variances was diagnosed via the plots of residuals vs. fitted values and vs. each predictor (see Figures S3 and S4). The residuals vs. fitted plot showed a curvilinear trend, suggesting the heterogeneity of the variance. To deal with this violation, the nlme package provides a standard class of variance function structures for specifying within-group variance models, e.g., fixed weights of a variance covariate (varFixed), constant variance (varIdent), exponential of a variance covariate (varExp). Since varIdent and varFixed were not applicable for a non-linear relationship between residual variance and covariates, varExp function was used as the variance was multiplied by an exponential function of the variance covariate Depth and an unknown parameter δ (Figures S8 and S9) [12].

To build a simple model, backward elimination strategy was applied [33]. Particularly, at first, maximum numbers of fixed effect variables were added in the model, i.e., COD, TKN, TP, pH, chlorophyll a, solar radiation, the log function of depth, the quadratic function of daytime and the interaction between daytime and solar radiation. After that, likelihood ratio tests were employed to test hypotheses about the fixed-effect variables in the LMM based on maximum likelihood estimation [33]. From that, non-significant predictors were identified and removed, i.e., COD (p-value = 0.173), TKN (p-value = 0.138), and TP (p-value = 0.495). The remaining variables, i.e., pH (p-value < 0.0001), chlorophyll a (p-value < 0.0001), solar radiation (p-value = 0.015), depth (p-value = 0.037), daytime (p-value < 0.0001), were kept in the final model as follows (Equation (5)).

$$\begin{array}{l}{\mathrm{LogDO}}_{\mathrm{ijk}} = -4.03 + {0.27 \times \mathrm{pH}}_{\mathrm{ijk}} + {0.0009 \times \mathrm{Chl}}_{\mathrm{ijk}} + 0.007 \times {\text{Solar.rad}}_{\mathrm{ijk}} - 0.59{ \times \mathrm{LogDepth}}_{\mathrm{ijk}} \\ + {0.51 \times \mathrm{Daytime}}_{\mathrm{ijk}} - 0.02{ \times \mathrm{Daytime}}_{\mathrm{ijk}}^2 - 0.001 {\times \mathrm{Daytime}}_{\mathrm{ijk}} \times {\text{Solar.rad}}_{\mathrm{ijk}} + 0.0001{ \times \mathrm{Daytime}}_{\mathrm{ijk}}^2\\ \times {\text{Solar.rad}}_{\mathrm{ijk}}{ + \mathrm{a}}_{\mathrm{k}}{ + \mathrm{a}}_{\mathrm{jk}}{ + \mathsf{\epsilon }}_{\mathrm{ijk}}\end{array}$$

(5)

where:

${\mathrm{a}}_{\mathrm{k}}$ ~ N(0,${\sigma }_{sc}^2$) with ${\widehat{\sigma }}_{sc}$ = 1.85 × 10⁻⁵,

${\mathrm{a}}_{\mathrm{jk}}$ ~ N(0,${\sigma }_p^2$) with ${\widehat{\sigma }}_p$ = 0.135,

${\mathsf{\epsilon }}_{\mathrm{ijk}}$ ~ N(0,${\sigma }^2$) with $\widehat{\sigma }$ = 0.116.

Surprisingly, none of the water-pollutant variables were in the final model, while pH and chlorophyll a appeared as significant predictors. The spatiotemporal effects were proved as daytime and depth were remained in the final model with the interaction between daytime and solar radiation. Regarding the random effects, random intercept for sampling campaign k, a_k, is normally distributed with mean 0 and very small variance (1.85 × 10⁻⁵)² while the random intercept for pond j within sampling campaign k, a_jk is normally distributed with mean 0 and much higher variance 0.135². From these variances, intraclass correlation coefficients were calculated, resulting in ICC_sc of 1.08 × 10⁻⁸ and ICC_pond of 0.58. Concerning the goodness-of-fit of the model, we obtained very high conditionals R² of 0.94, which suggested both fixed and random effect variables providing considerable potential in predicting the oxygen level within the pond system.

3.3. Model Evaluation

To assess the predictive performance of the model, a scatter plot of observed vs. predicted values was drawn (see Figure 5). Testing of the model predictions against observed data demonstrated that the model was capable of capturing the variation of oxygen level within different ponds. Indeed, except for the lowest concentrations of oxygen being slightly underestimated, the model predictions deviated less than an order of magnitude, inferring high model accuracy and consistency.

Figure 6 showed the MAE values between predicted and observed oxygen concentrations at two different depths in each individual pond and the whole system. Generally, relatively low MAEs indicated the fairly high accuracy of the mixed model. Apart from the high MAEs of FP1, this measure of error was lower than 1.3 mg O₂·L⁻¹ in the other ponds. In fact, 26 out of 30 highest values of MAE were from the samples collected at FP1. The other four values belonged to the samples at MP2, causing this pond having the second highest values of MAE. It is noteworthy that both of these ponds had the highest fluctuation of oxygen level at the two depths.

3.4. Diurnal Dissolved Oxygen Profile

The hourly variation of oxygen concentration was illustrated in Figure 7. Particularly, the graph depicted the mean observed and predicted oxygen concentrations with their equivalent values of MAE. These values were averaged over four ponds at two different depths. At the surface, the observed oxygen concentration gradually increased from 9:00 a.m. onwards and reached its peak, around 9 mg O₂·L⁻¹, from noon to 5:00 p.m. On the other hand, the increasing trend started later at the bottom from around 12:00 p.m. and then the oxygen level surprisingly boosted to the same concentration at the surface layer during the last two hours. In fact, during this time of the first sampling campaign in FP1, a heavy rain was recorded, causing abnormal oxygen concentrations at the bottom layers, up to 16 mg O₂·L⁻¹. Comparing to the collected data, the mixed model was able to describe relatively precisely the diurnal oxygen variation, especially at the surface where the average MAE was only 0.80 mg O₂·L⁻¹. Similarly, except for the last two hours when the oxygen level was greatly underestimated, the average MAE values for the prediction of the samples at the bottom layers were very low, i.e., 0.52 mg O₂·L⁻¹.

3.5. Vertical Dissolved Oxygen Profile

Figure 8 demonstrated the vertical dissolved oxygen profile of the pond system, in which predicted and observed oxygen concentrations with the corresponding values of MAE were averaged over four ponds at different depths. As showed in the figure, the mixed model predicted the vertical pattern of oxygen reasonably well as it captured the sudden drop of oxygen level from epilimnion to hypolimnion layer, which can be a consequence of algal stratification due to light limitation [5]. There were three abnormally high oxygen concentrations at depths 65, 85, and 115 cm, which were recorded at FP1 in the heavy-rain conditions, causing high MAEs of around 3 mg O₂·L⁻¹.

4. Discussion

4.1. Spatiotemporal Influences on the Oxygen Dynamic

Spatiotemporal autocorrelations can be seen as both an opportunity and a challenge as it can provide useful information for inference of process from their patterns [22]. In fact, to account for these autocorrelations, we implemented the first application of mixed model in WSPs, which underlined the variations of oxygen concentrations with respect to time, and space in a high-altitude WSP. The LMM (Equation (5)) showed the negative impact of depth (regression coefficient = −0.59) and the diurnal effects with the shape of a downward opening parabola as the coefficient of daytime² is −0.02. These results agreed with the observed data from the three sampling campaigns as shown in Figure 7 and Figure 8 as well as the high value of conditional R² of 0.94. Interestingly, none of the water-pollutant variables, i.e., COD, TKN, and TP, was included as the fixed-effect variables in the final model. The reason could be because of high oxygen concentrations as a result of the enhanced algal photosynthesis, which was accelerated by strong solar radiation in this meridional high-altitude WSP system. Compared to this high production of DO, the amount of oxygen, which was consumed for bacterial mineralization and nitrification process, appeared to be inconsiderable. Indeed, the remaining fixed-effect variables are related to only the algal photosynthesis and their positive coefficients, i.e., pH (0.27), chlorophyll a (0.0009), and solar radiation (0.007), also support the substantial influence of algal photosynthesis on the spatiotemporal variations of oxygen concentration.

Moreover, the spatiotemporal variations of DO were highlighted via the random effects. In particular, the very low variance of the random sampling-campaign effects (${\sigma }_{sc}^2$) leading to almost zero ICC_sc suggested that there was nearly no dissimilarity in the oxygen variation among the three sampling campaigns. The reason of this similarity can be explained by the fact that the sampling campaigns were conducted within one dry season of Ecuador when a very widely fluctuated oxygen level can be observed over the course of one day as showed in Figure 7. In Cuenca, daylight lasted from 6:00 a.m. to 6:00 p.m. with peak of solar radiation between 1:00 p.m. and 3:00 p.m. promoting the maximum rate of algal photosynthesis [35]. During this peak period, an extremely high oxygen level was recorded, i.e., more than 20 mg O₂·L⁻¹ in our sampling campaigns and up to 39 mg O₂·L⁻¹ in the sampling campaign of Alvarado [35] in this pond system. These abnormally high oxygen levels can be induced by vast light intensity at high altitude, up to 1500 W·m⁻², and high algal biomass above 420 µg Chl a·L⁻¹ near the surface of the FPs. This high algal biomass can be a result of intensive solar radiation at high altitude. In fact, in a WSP located at an altitude of 2675 m in Mexico, Pearson et al. [36] also found extremely high levels of chlorophyll a, up to 1500 µg Chl a·L⁻¹. This algal overgrowth can generate a supersaturated DO condition during the day, but, on the other hand, depletes the oxygen level due to their respiration during the night [37]. As such, a vast fluctuation of the pond performance can be found between early morning and mid-afternoon in a high-altitude WSP.

Contrast to the very low ICC_sc, the relatively high value of ICC_pond of 0.58 suggests that both the variance of the random effects associated among the ponds (${\sigma }_p^2$) and the variation of the oxygen concentration within a pond (${\sigma }^2$) were considerable. Indeed, the variations associated among the ponds derived from the difference in pond performance across the two flow lines. More specifically, FP1 received around 20% higher pollutant loadings than FP2, especially organic matter. In fact, their average surface organic loadings were up to 250 and 185 kg·ha⁻¹·day⁻¹, respectively, while the recommended limitations were 240 kg·ha⁻¹·day⁻¹ for WSPs at tropical and subtropical regions and only 200 kg·ha⁻¹·day⁻¹ for WSPs at altitudes above 2400 m a.s.l. [13,38]. These high loadings were associated with the sludge accumulation, which was also the reason of unavailable data at the bottom of location 1 and 2 in FPs. According to the study of Alvarado et al. [39] in this pond system, the sludge volume of the FPs reached up to 34% of the pond volume, which substantially reduced its active volume.

Regarding the variation within a pond, in contrast to a relatively homogeneity of oxygen level in the MPs, higher concentrations were found at the central area of the FPs. This difference can be caused by the fact that the central area of the FPs had less sludge accumulation and flow turbulence, which promoted high density of algae located at this region [39]. On the other hand, the bathymetries in the study of Alvarado et al. [16] on these MPs showed only a slight sludge layer growth in the maturation ponds that can be assumed negligible regarding the pond hydraulics.

4.2. Model Evaluation

The mixed model proved its ability to capture the very dynamic variation of oxygen level in a high-altitude pond system. This is supported by very high goodness-of-fit, a fair agreement between predicted and observed values (Figure 7 and Figure 8), and low values of MAEs, which were normally smaller than 10% of the predicted data, except for the abnormally high DO observed in FP1. This abnormality was generated due to the heavy rain with high wind speed, above 5 m·s⁻¹, leading to high turbulent flow that homogenized the water column, creating very high concentrations of DO, up to 16 mg O₂·L⁻¹. These values caused poor predictive performance of the mixed model during the last two hours of the first sampling campaign. This underestimation can be explained by the sensitivity of mixed models to abnormal observations [40] and the fact that the three sampling campaigns were conducted within one dry season of Ecuador. As the photosynthetic activity of algae changes in response to the seasonal changes in environmental conditions [24,41], there is a need for additional sampling campaigns in the rainy season. From that, the interpolation of this mixed model can be reliably performed within a larger range of observation.

Concerning the model applicability, a question should be raised, related to the usefulness of this mixed model in terms of predicting oxygen concentrations compared to previous mechanistic models. Firstly, when encountering with the complex interactions of multifaceted factors in ecological systems, a statistical model can be preferred over a mechanistic model due to its simplicity [42]. Especially noteworthy is that the WSP system as an open natural system should be considered as a complicated assemblage of different processes and inputs, hence, its mechanistic models are nearly always overparameterized [43]. In fact, to simplify the models, some important processes were neglected in Kayombo, Mbwette, Mayo, Katima and Jorgensen [6] and Banks, Koloskov, Lock and Heaven [7], e.g., nutrient uptake of bacteria and algae, nitrification/denitrification, and air/water exchange. More importantly, numerous parameters applied in mechanistic models have been taken from different systems and model assumptions based on external characteristics might have to be applied. This approach of artificially assigning values to parameters can lead to biased results, which significantly deviates from real outputs [38]. Indeed, several simplifying assumptions along with geometric approximations caused the disagreements between computational fluid dynamics (CFD) model outputs and experimental results in the research of Alvarado, Sanchez, Durazno, Vesvikar and Nopens [39]. Finally, the large number of experimental data needed for the validation of mechanistic model is another major constraint.

4.3. Insights for Oxygen Regulation in WSPs

According to previous guidelines of Pearson et al. [44], and Mara and Pearson [45], DO was identified as an additional parameter for effluent quality monitoring and evaluation of pond performance, as it is determined not necessary to monitor this extra parameter for routine monitoring and evaluation. Nevertheless, this is totally not the case in conventional biological WWTPs using activated sludge. Oxygen levels are considered a key parameter in the operation of such a plant and DO control is of primary importance in activated sludge processes [46]. Particularly, DO concentration should be sufficient for aerobic microorganisms to degrade organic matters and convert ammonium to nitrate, yet it is not excessive to deteriorate the sludge formation, which can lead to the problem of sludge settleability. In addition, as air supply accounts for the largest portion of the operational budget of these plants, proper DO control can save enormous energy costs compared to uncontrolled systems [47]. In fact, to optimize the operation of WWTPs, DO control via aeration is one of a very limited number of manipulatable variables so that has been the subject of extensive research since the 1970s [48]. In contrast to this maturity, the omission of DO control in WSPs can be a neglected reason causing many common problems in pond operations, such as organic overloading, odor problems, and algal overgrowth. This lack of proper control and poor operation was reported as one of the main reasons for the under-performance of eight high-altitude WSP systems in Mexico [49]. Given these points, pond engineers need a more advanced strategy for O&M in general and DO control in particular. To do so in conventional WWTPs, instrumentation, control, and automation (ICA) as an advanced tool for system control has been long developed and is now well-recognized as an integrated part of plant operation while this technology has been very new in pond operation up to now [48]. This limitation is because, for small-scale ponds, the O&M tends to be neglected due to financial reasons while, in large-scale systems, pond engineers are still more comfortable with the traditional procedures [5]. However, since effluent discharge standards become increasingly stricter and levy charges for plant performance failures can sharply increase the O&M costs, advanced control appears to be a more economic and reliable method. As such, we propose following practical recommendations on DO monitoring and regulation in WSPs, based on ICA technology and the findings from our research.

Firstly, the flow behavior of a pond should be taken into account in the DO control as it can considerably affect oxygen variations. While homogenous oxygen levels are normally observed in completely mixed systems, the DO profile in plug-flow-like ponds, which reflects oxygen uptake rate (OUR), can vary greatly along the systems [50]. Different from activated sludge systems where these two ideal regimes are commonly applied, the mixing behavior in WSPs is more complicated as being affected by many factors, such as inlet/outlet configurations, wind, sludge accumulation, and baffles [51]. As such, in order to evaluate DO profile and flow performance, we suggest that the locations in the grid scheme in Figure 1, which cover the three essential areas of a pond, i.e., influent, middle, and effluent, are needed to be sampled. From that, we can identify the flow stratification and recognize dead spots, stagnant areas, and possible sludge accumulation as it was the influent area in our case. It is noteworthy that baffled ponds, where the hydraulic regime is similar to the ideal plug flow, can offer more flexibility for DO control but also greater challenges [52]. In plug-flow systems, different DO set-points can be chosen along the reactors, providing independent zones that can be used for different purposes, such as aerobic zones for nitrification and anoxic zones for denitrification. Such a configuration leads to a more complex control system, which needs to include ammonia and nitrate sensors like in Kallby WWTP (Sweden), which has ten zones with a pre-denitrification configuration [47]. This potentiality in baffled ponds to increase the nutrient removal can be exploited via optimal control and design to setup the most appropriate measurements and their corresponding set-points. Such a study is surely guaranteed, but still needs to be done.

Secondly, pond engineers also need to evaluate the vertical DO profile which is a result of algal stratification reflecting light availability along the water depth. Generally, light and photosynthetic activity can extend down to the bottom of shallow MPs where clean and less turbid water is located, while, with higher concentration of suspended solids, this extension for FPs is only 20 to 30 cm [5]. This light attenuation can be very important for pond performance. The location of the sudden drop of oxygen defines the volume of anoxic and anaerobic area in FPs (see Figure 8). Vigorous mixing can carry oxygenated water from upper aerobic layers to the bottom area, which limits the extent of methane fermentation, leading to acidic conditions and odor release [53]. This common problem of shallow FPs occurred during the last two hours in the first sampling campaign when the heavy rain and strong wind homogenized the oxygen level within the water column in FP1 (see Figure 7). The occurrence of such a disturbance is the major incentive of system control. Naturally, WSPs are relatively resilient to disturbances as a result of large volume and thus high HRT. However, besides encountering a large variation of wastewater influent regarding both its composition and flow rates, WSPs, especially the systems at high altitude, also have to deal with the hourly, daily and seasonal changes in climatic conditions. A traditional way of tackling this issue was to build a larger volume as it was a suggestion of Juanico et al. [13] for high-altitude WSPs. Compared to ICA technology, it is not the best economic solution, as overly conservative designs can inflate capital and O&M costs of the plants and is not an optimal choice for pond upgrading [52]. Regarding MPs, two key mechanisms of disinfection process involve photo-oxidation, which essentially relies on the presence of DO and pH [54]. Indeed, it was concluded in Curtis et al. [55] and Dixo et al. [56] that sunlight-mediated disinfection can only have a considerable impact on fecal coliforms in the case of high DO and pH present. Interestingly, high algal biomass generates high DO and pH, and reduces the light penetration in the ponds since algae contain large amounts of pigments that can block sunlight [57]. Hence, in MPs with low concentrations of other light absorbers, such as gilvin and tripton, the information from the DO profile showing the algal distribution and light penetration can facilitate the evaluation of pond disinfection efficiency.

Moreover, as 80% of the oxygen source in ponds originates from algal photosynthesis, pond engineers should keep in mind its diurnal variation. The period of daylight should be recorded, as solar radiation is the main energy provider for the ponds. It is also recommended that the sample of DO should be collected periodically at a minimum of three moments in a day, i.e., sunrise, noon, and sunset. Extra care should be given in high-altitude pond systems since high light intensity promoting algal overgrowth can generate supersaturation of oxygen during daylight but drain out oxygen when the light diminishes. This can lead to overload and violate the discharge permit; hence, extra aeration may be needed during the night. Furthermore, as not only algal photosynthetic activities but also other characteristics of ponds, such as nutrients, bacterial levels, and dissolved organic matter changes seasonally [24,25], it is advised that the oxygen profile should be recorded and compared between different seasons during a year for an accurate depiction of pond performance.

5. Conclusions

The first application of a linear mixed effect model highlighted the spatiotemporal variations of oxygen level in WSPs, which were enhanced by severe meteorological conditions at high altitude in this study. Particularly, a substantial diurnal variation was observed from zero to above 20 mg O₂·L⁻¹, which can be a result of algal overgrowth as it was expedited by the intensive solar radiation at 2400 m a.s.l. This algal bloom generated supersaturation of oxygen during the day but drained out oxygen via their respiration during the night. The critical role of algae in the oxygen temporal dynamics was also emphasized in the final model, as all the remaining fixed-effect variables were associated with the algal photosynthesis. Despite being designed in the two parallel flow lines, different ponds exhibited different spatial patterns of oxygen dynamics as a result of numerous factors, such as flow behavior, sludge accumulation, algal distribution, influent fluctuation, and pond function. In the mixed model, this spatial variation was indicated via the high variance of the random effects associated among the ponds, ICC_pond of 0.58.

From these findings, together with the fact that DO control in WSPs is overlooked in the previous guidelines of pond O&M, some practical recommendations are given. Particularly, hydraulic performance should be taken into account in DO control, which can be very advantageous for baffled ponds to optimize nutrient removals by optimal control and design to setup proper measurements and their corresponding set-points. Pond operators should also pay more attention to the vertical DO profile, which reflects algal distribution and light penetration. As these factors play an important role in pond functions, the information from the vertical DO profile can facilitate the evaluation of pond performance. Especially noteworthy in the case of high-altitude WSPs is that the variation of climatic conditions should be recorded, i.e., light intensity, cloudiness, precipitation and air temperature. Unusual disturbances from extreme climate can lead to high levy costs for discharge violations, which has been proved from conventional WWTPs that can be mitigated by the application of advanced system control, i.e., ICA technology. More importantly, since the mixed model proved its ability to cope with high autocorrelations of data in pond systems, and from that provided more useful information on spatiotemporal patterns, we recommend this advanced statistical technique as an effective tool for better understanding and simulation to pond engineers and researchers.