OPTIMIZING COLD-CLIMATE TREATMENT WETLAND PERFORMANCE: INTEGRATING MACHINE LEARNING ANALYSIS WITH EXPERIMENTAL EVALUATION OF CAREX UTRICULATA'S FUNCTIONAL RESILIENCE by Erik Joseph Groenenberg B.A.Sc., University of British Columbia, 2017 THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE IN ENGINEERING UNIVERSITY OF NORTHERN BRITISH COLUMBIA April 2025 © Erik Joseph Groenenberg, 2025 Abstract Treatment wetlands experience significant seasonal temperature variations that affect biological treatment processes through complex interactions between plants, microorganisms, and environmental conditions. This project employs a dual methodological approach to optimize cold-climate treatment wetland design and operation. First, we developed interpretable machine learning models using data from 27 published studies and 118 treatment wetlands to predict effluent temperature (RMSE = 0.7°C) and ammonia (RMSE 2.9 – 9.3 mg/L) and organic matter (RMSE 11 – 44 mg/L) concentrations in different wetland configurations. Model interpretation revealed that influent temperature is the dominant predictor of effluent temperature followed by air temperature, with hydraulic loading rate modifying this relationship. For contaminant removal, we found that plant species' coldtemperature benefits are most pronounced in saturated systems, while their impact is marginal in unsaturated configurations. Building on these insights, we conducted controlled microcosm experiments with Carex utriculata in batch-operated wetlands across temperature conditions (23°C and 5°C) and after harvesting. C. utriculata maintained superior organic matter removal in cold conditions (86±5% versus 73±9% for unplanted systems) and continued to outperform unplanted systems even after harvesting. Evapotranspiration-driven water level reduction in warm conditions improved porosity maintenance by 2-3% compared to unplanted systems. Our integrated approach demonstrates that C. utriculata provides resilient cold-temperature treatment benefits while offering operational advantages through reduced clogging potential and maintained performance post-harvest, addressing key challenges identified in machine learning analysis of published literature. ii Table of Contents Abstract ................................................................................................................................... ii Table of Contents ................................................................................................................... iii List of Tables ......................................................................................................................... vii List of Figures ...................................................................................................................... viii Glossary ................................................................................................................................ xiii Acknowledgements ...............................................................................................................xiv Chapter 1 Introduction.................................................................................................... 15 1.1 Chapter overview ................................................................................................... 15 1.2 Background ............................................................................................................ 15 1.3 Problem statements ............................................................................................... 19 1.3.1 Limitations of temperature inclusion and model interpretation of treatment wetlands ............................................................................................................................ 19 1.3.2 Batch operation of C. utriculata treatment wetlands through seasonal changes . 19 1.4 Objectives ............................................................................................................... 20 1.5 Thesis organization ................................................................................................ 21 1.6 Scientific outcomes and outputs ........................................................................... 21 Chapter 2 Interpretable Machine Learning Reveals Cold Climate Dynamics of Treatment Wetlands .............................................................................................................. 23 2.1 Introduction............................................................................................................ 23 iii 2.2 Methods .................................................................................................................. 25 2.2.1 Data collection ..................................................................................................... 25 2.2.2 Model construction .............................................................................................. 27 2.2.3 Interpretation........................................................................................................ 28 2.3 Temperature results and discussion ..................................................................... 30 2.3.1 Data collection ..................................................................................................... 30 2.3.2 Model performance .............................................................................................. 31 2.3.3 Model interpretation ............................................................................................ 32 2.4 Treatment results and discussion ......................................................................... 37 2.4.1 Data collection ..................................................................................................... 37 2.4.2 Preliminary modelling and data division ............................................................. 39 2.4.3 Model performance .............................................................................................. 42 2.4.4 Model interpretation ............................................................................................ 42 2.5 Limitations ............................................................................................................. 55 2.6 Conclusion .............................................................................................................. 57 Chapter 3 Using Machine Learning Insights to Target Research on the Seasonal Effects of Plants in Batch Treatment Wetlands .................................................................. 59 3.1 Introduction............................................................................................................ 59 3.2 Insights from machine learning on the effects of plants .................................... 59 3.3 Research questions driven by machine learning insights .................................. 61 iv 3.4 Methodological decisions ...................................................................................... 63 3.5 Expected outcomes, connections, and significance ............................................. 65 Chapter 4 C. utriculata Improves Treatment in Three-Day Batch Operated Treatment Wetlands in Warm, Cold, And Harvested Conditions ..................................... 66 4.1 Introduction............................................................................................................ 66 4.2 Materials and Methods.......................................................................................... 68 4.2.1 Experimental setup .............................................................................................. 68 4.2.2 Synthetic wastewater composition ...................................................................... 69 4.2.3 Operation, sampling, and measurements ............................................................. 70 4.3 Results and discussion ........................................................................................... 72 4.3.1 Evapotranspiration ............................................................................................... 72 4.3.2 Chemical Oxygen Demand .................................................................................. 73 4.3.3 Ammonia ............................................................................................................. 77 4.3.4 Phosphorus........................................................................................................... 79 4.3.5 Porosity ................................................................................................................ 80 4.4 Conclusion .............................................................................................................. 82 Chapter 5 Conclusion ...................................................................................................... 83 5.1 Background ............................................................................................................ 83 5.2 Joint results of machine learning and microcosm experiments ........................ 83 5.3 Cold climate design and operation ....................................................................... 86 v 5.4 Limitations and recommendations for future work ........................................... 88 References ............................................................................................................................... 92 Appendix A : Modelling datasets ................................................................................. 119 Appendix B : Supplementary model interpretation plots ......................................... 120 vi List of Tables Table 2-1: P-K-C* variables which correspond to the mean treatment wetland features in our machine learning data sets. These are used to provide predictions of effluent concentrations to compare with model interpretation plots......................................... 30 Table 2-2: Features considered for training the temperature model and the reasons for including or not including them................................................................................... 31 Table 2-3: Features considered for inclusion in the treatment model training data, including whether they were included and the rationale. ............................................................ 38 Table 2-4: Mean values of features correlated with saturation type in the treatment data set. as well as how many unique treatment wetlands were collected within each saturation type. ............................................................................................................................. 40 Table 2-5: The range, mean, and standard deviations of numeric features in the separated treatment datasets......................................................................................................... 41 Table 2-6: Fit metrics for each model and algorithm that is predicting effluent concentrations. ..................................................................................................................................... 42 Table 4-1: Synthetic wastewater composition ......................................................................... 70 Table 4-2: Air and water temperature throughout the experiment, depending on the location of the microcosms. Values are listed as mean ± standard deviation (min – max) to show variability. .......................................................................................................... 72 Table A-1: Papers used in temperature dataset. ..................................................................... 119 Table A-2: Papers used in treatment datasets ........................................................................ 119 vii List of Figures Figure 2-1: Effluent temperature predictions plotted versus actual effluent temperature for each model algorithm tested, as well as their R2 and RSME goodness of fit scores. .. 32 Figure 2-2: Ranked importance of features and interactions between features for the prediction of effluent temperature. SHAP values show the overall importance of each feature. Interaction values show how much variance in predictions is attributed to interactions of either a particular feature (interaction strength) or a specific interaction (shown with “:” between the interacting features). ..................................................... 33 Figure 2-3: Accumulated local effect of each feature in the model. The marks along the x-axis indicate data point presence, giving an approximation of data distribution. ............... 34 Figure 2-4: The accumulated local effect of delta on the predicted effluent temperature depending on the specific algorithm (right), and the accumulated local interaction effect between delta and HLR for the Cubist (top-left) and XGBoost (bottom-left) algorithms. The colour legend for the interaction effects is shown in the top-right. Presence of training data points are shown as black dots as well as black lines along the x and y axis. ........................................................................................................... 35 Figure 2-5: ALE of delta on the predicted effluent temperature for different algorithms (separate plots) and saturation types (different coloured lines). Training data presence is noted by marks along the x and y axis, coloured by saturation type. ...................... 36 Figure 2-6: SHAP values show the relative importance of each feature depending on if ammonia or organic matter is being predicted (top labels, first line) and if the training data is from saturated or unsaturated treatment wetlands (top labels, second line). viii Coloured dots show SHAP values of specific algorithms, while the box plots show their median importance and distribution. ................................................................... 43 Figure 2-7: PDP and ICE plots show how organic matter influent (x-axis) [mg/L] changes the predicted organic matter or ammonia effluent (y-axis) depending on the algorithm (colored lines) and saturation type (separate plots). .................................................... 44 Figure 2-8: PDP and ICE lines show how ammonia influent [mg/L] influences the predicted ammonia effluent for saturated (left) and unsaturated (right) treatment wetlands depending on the algorithm (coloured line). Unsaturated vertical wetlands had no training data with influent ammonia greater than 100 mg/L. ...................................... 45 Figure 2-9: PDP and ICE plots of how HLR [L/m2/d] influences each algorithms predicted organic matter or ammonia effluent, depending on the saturation of the treatment wetland......................................................................................................................... 46 Figure 2-10: PDP and ICE plots showing how the minimum media diameter [mm] influences the predicted ammonia and organic matter concentrations for each algorithm, depending on the saturation of the treatment wetland. ................................................ 47 Figure 2-11: PDP and ICE plots showing how the depth of the treatment wetland bed [m] influences the predicted ammonia and organic matter concentrations for each algorithm, depending on the saturation of the treatment wetland. .............................. 48 Figure 2-12: PDP and ICE plots showing how effluent temperature [°C] influences the predicted organic matter and ammonia effluent concentrations for each algorithm, depending on the saturation of the treatment wetland. ................................................ 49 ix Figure 2-13: PDP plots showing how effluent temperature [°C] influences the predicted effluent organic matter concentrations for all saturated systems (continuous and batch), depending on the algorithm (from left to right), and plant status (coloured lines). ........................................................................................................................... 50 Figure 2-14: PDP plots showing how effluent temperature [°C] influences the predicted effluent organic matter concentrations for continuous flow saturated treatment wetlands (no batch), depending on the algorithm (from left to right), and plant status (coloured lines). ........................................................................................................... 51 Figure 2-15: PDP plots showing how effluent temperature [°C] influences the predicted effluent organic matter concentrations for vertical unsaturated systems, depending on the algorithm (from left to right), and plant status (coloured lines). ........................... 52 Figure 2-16: PDP and ICE plots showing how effluent temperature [°C] influences the predicted effluent ammonia-N concentration, depending on the saturation (upper and lower plots), algorithm (from left to right), and plant status (coloured lines). ............ 53 Figure 2-17: Plot of influent (yellow) and effluent (blue) water temperatures versus air temperature for all outdoor buried systems in the temperature dataset. ...................... 54 Figure 2-18: PDP and ICE plots showing how the treatment wetland age [years] affects the predicted organic matter and ammonia-N effluent concentrations, depending on the saturation and algorithm. ............................................................................................. 55 Figure 4-1: From left to right: microcosms in the EFL (batches 0 – 29), the cold room (batches 29 – 50), and after being harvested in the cold room (batches 50 – 70) ....... 71 x Figure 4-2: Evapotranspiration (for planted units) and evaporation (for unplanted units) throughout each batch of the experimental period. The first dotted line marks the start of cold conditions, while the second marks when plants were harvested. .................. 73 Figure 4-3: NH3-N, COD, and PO4-P effluent concentrations (left) and mass (right) for planted (red) and unplanted (blue) microcosms. The first line marks the start of cold room conditions, and the second line marks the time at which C. utriculata was harvested. ..................................................................................................................... 74 Figure 4-4: pH in the effluent of each batch. The first line marks the start of cold room conditions, and the second line marks the time at which C. utriculata was harvested. ..................................................................................................................................... 75 Figure 4-5: Effective porosity of each batch, based on the volume to fill the drained microcosms to a set water level. The first line marks the start of cold room conditions, and the second line marks the time at which C. utriculata was harvested. ................. 81 Figure 5-1: PDP curves of the effect of plants dependant on temperature AFTER Chapter 4 microcosm data was added. Results are in between the original saturated model and the continuous flow only model. Harvesting shows a larger effect than previously predicted. ..................................................................................................................... 84 Figure 5-2: PDP curve for the effect of temperature depending on the whether the treatment wetland is fill-and-drain (batch) operated or permanently saturated. Model shows results AFTER microcosm data added. ....................................................................... 85 Figure 5-3: PDP curves AFTER microcosm data added for the effect of temperature, depending on the organic matter parameter for saturated systems. COD is more xi consistently greater than sCOD in the revised plots, especially at colder temperatures. ..................................................................................................................................... 86 Figure B-1: Pearson correlation and data distribution for features in temperature dataset. X and Y axis labels are shown along the central diagonal. ........................................... 120 Figure B-2: PDP and ICE curves for the effect of loading interval [hours] on ammonia and organic matter concentrations .................................................................................... 120 Figure B-3: PDP curves for the effect of temperature depending on the organic matter parameter for saturated systems................................................................................. 121 Figure B-4: PDP curve for the effect of temperature depending on the whether the treatment wetland is fill-and-drain (batch) operated or permanently saturated. ........................ 121 Figure B-5: PDP curve for the effect of temperature depending on the whether the treatment wetland is fill-and-drain (batch) operated or permanently saturated. AFTER microcosm data added. .............................................................................................. 121 Figure B-6: PDP and ICE points for the effect of plants on ammonia and organic matter concentrations ............................................................................................................ 122 Figure B-7: PDP and ICE curves for the effect of recirculation rate [fraction of inflow] on ammonia and organic matter concentrations ............................................................. 122 xii Glossary ALE accumulated local effects BOD5 five-day biochemical oxygen demand C* background concentration CBOD carbonaceous biochemical oxygen demand COD chemical oxygen demand EFL enhanced forestry laboratory GBM gradient boosting machine HDPE high density polyethylene HLR hydraulic loading rate HRT hydraulic retention time ICE individual conditional expectation k areal removal rate constant NH3 ammonia Org-N organic nitrogen P tanks-in-series constant PAR photosynthetically active radiation PDP partial dependance plot PO4 phosphate RF random forest RSME roost square mean error sCOD soluble chemical oxygen demand SHAP Shapley additive explanations SVM support vector machines tCOD total chemical oxygen demand TN total nitrogen XGBoost extreme gradient boosting xiii Acknowledgements This thesis would not exist without the assistance of many people and organizations. To my supervisors Dr. June Garcia-Becerra and Dr. Oliver Iorhemen, I thank you for your excellent advice, patience, and encouragement over my time at UNBC. To my committee members Dr. Cindy Smith and Dr. Theresa Adesanya, thank you for your input, excitement, and willingness to help. Many hours were spent in the Enhanced Forestry Lab, where Doug Thompson and Kennedy Boateng were extremely generous with their assistance, EFL space, and time to chat about the project and life. Thanks to Fryderyk Packowski of the Hydrology Lab and Spark Labs as well as Caleb for lending and supporting our use of probes for the project, and my father-inlaw Andy for generously building and providing a probe monitoring system. Thank you to all the students in both June and Oliver’s research groups for their support, especially Mario, Bridget, Raul, Carla, and Sofia for their assistance, advice, and encouragement. A massive thank you to R. Radloff and Associates Inc., who encouraged me to pursue a master’s degree, allowed me to take two years of leave from work, and provided financial support. Thanks also to UNBC: the generous BC Graduate Scholarship, the UNBC Research Project Award, and the Green Travel Grant. Thanks also to the many people I met at the International Water Association who expressed great interest in my work, offered data, and welcomed me into the scientific community. Above all I thank my wife Johanna for her constant support and encouragement throughout. xiv Chapter 1 Introduction 1.1 Chapter overview Constructed treatment wetlands represent a broad range of engineered and purpose-built systems for the treatment of wastewater, inspired by natural wetlands. Treatment wetlands pass wastewater through (subsurface) or above (free-water) granular or soil media planted with wetland species, where biological, chemical, and physical processes treat the water (Kadlec & Wallace, 2009). Unsurprisingly this results in a myriad of design and operational possibilities, which only increase in complexity once variations in wastewater composure, climate, and natural plant growth are considered. This thesis focuses on subsurface treatment wetlands using machine learning models trained on data published by others (Chapter 2) and experimental microcosms (Chapter 4) to investigate the dynamics between design and operation, temperature changes, plant effects, and effluent quality. 1.2 Background As treatment wetlands are outdoor systems taking advantage of natural processes, they can be affected by seasonal changes. This thesis focuses on cold climates, defined by the KöppenGeiger climate classification system as having a mean air temperature below 0°C for at least one month of the year (Cui et al., 2021). In these climates plant dormancy and senescence become relevant, and the wetland bed reaches colder temperatures that can lower treatment efficiency. The impact of water temperature depends on the design and operational parameters which determine the environmental conditions within the wetland. Water temperature changes in the wetland can impact microbial activity, which is the primary means of treatment (Samsó & 15 Garcia, 2013). Colder temperatures reduce microbial affinity for substrates (Nedwell, 1999), however increasing bioavailable reaction-limiting components can compensate for temperature decrease (C. R. Allen et al., 2023; Rusten et al., 1995; Salvetti et al., 2006). Treatment wetland designs have very different oxygenation abilities (Nivala et al., 2013), altering their temperature response. Water temperature itself is not directly correlated with seasons or climate but controlled by the energy balance (Equation 1) (Kadlec & Wallace, 2009). Influent temperature has shown to have the greatest influence (Grebenshchykova et al., 2023). While influent temperature has a linear relationship with air temperature in warmer conditions (>10°C) (Mietto et al., 2015) this can diverge as air temperatures lower. Outside design elements such as wastewater collection system length and snowmelt or groundwater infiltration and inflow (Panasiuk et al., 2022; Prost-Boucle et al., 2015a; Rusten & Ødegaard, 2023; Vidal et al., 2023), can all alter the influent temperature, requiring it to be considered separately from air temperature or season. Further modifying the wetland temperature are design parameters such as the hydraulic loading rate (Equation 1), air inflow rate of different designs and operations, and insulative cover (Agnieszka & Anna, 2013; Grebenshchykova et al., 2023; Prost-Boucle et al., 2015a). 𝑇𝑇𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 = �𝑅𝑅𝑁𝑁 + ℎ𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐 (𝑇𝑇𝑎𝑎𝑎𝑎𝑎𝑎 − 𝑇𝑇𝑡𝑡𝑡𝑡𝑡𝑡 𝑜𝑜𝑜𝑜 𝑏𝑏𝑏𝑏𝑏𝑏 ) + 𝑇𝑇𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝜌𝜌 𝑐𝑐𝑝𝑝 𝐻𝐻𝐻𝐻𝐻𝐻 − 𝜆𝜆𝑚𝑚 𝜌𝜌 𝐸𝐸𝐸𝐸 − 𝐺𝐺� �𝜌𝜌 𝑐𝑐𝑝𝑝 𝐻𝐻𝐻𝐻𝐻𝐻� Equation 1: Effluent temperature as a function of energy balance RN: Net radiation absorbed by the wetland hconv: Convective heat transfer coefficient [MJ/m2/d/°C] λm: Latent heat of vaporization of water 16 ρ: Density of water ET: Water lost to evapotranspiration G: Net loss to ground HLR: Hydraulic loading rate Seasonal changes of air temperature more directly affect the plants in the wetland. Nutrient uptake (Thorén et al., 2004a), root oxidation (Taylor et al., 2011), and evapotranspiration of plants (Beebe et al., 2014) all influence the effluent quality and vary with season. In summer, evapotranspiration draws up water from the treatment wetland, altering the hydraulics and in some cases drying up part of the treatment wetland if greater than the hydraulic loading rate (Beebe et al., 2014; Nivala et al., 2022). As colder temperatures slow evapotranspiration and the metabolism of microbial species, radial oxygen loss through plant roots has been found to enhance aerobic treatment efficiency, with Carex utriculata more effective than any other species at very low hydraulic loading rates in batch operated microcosms (Taylor et al., 2011). Contrasting with this benefit, there is a need to cut and remove the above ground plant stalks periodically, often during fall. The nutrients plants take up into their biomass provide a valuable treatment mechanism, but these nutrients return to the treatment wetland if allowed to decompose over winter (Vymazal, 2007). Wastewater composition, treatment wetland operation and design, and maturation over time also alter treatment processes. Hydrolysis of particulate organic matter is reduced in cold temperatures (Y. Huang et al., 2005), resulting in particulate matter build up and capture by plant roots (Karathanasis et al., 2003; Samsó & Garcia, 2013). Variable water levels brought about by evapotranspiration or fill-and-drain batch operation may alter biofilm exposure to 17 air in seasonally varied hydraulic conditions, altering root (J. Zhang et al., 2024) and biofilm (Hua et al., 2014a) formation. These all tie into longer term biofilm growth and particulate build up that alters the effective porosity and hydraulic conductivity of treatment wetlands (H. Wang et al., 2021). Studies have worked to integrate these many treatment wetland variables into cohesive models. Mechanistic approaches successfully simulate our understanding of the many processes within treatment wetlands (Langergraber & Šimůnek, 2012; Samsó & Garcia, 2013) but can be impractical to use for design purposes due to their complexity and many input variables. However, by calibrating these models on several real treatment wetlands, simulations have generated many useful general insights about the processes within wetlands and how water temperature and other design parameters jointly alter effluent quality (Martí et al., 2018; Pucher & Langergraber, 2019). As data gets more available, machine learning models have exploded in usage in biological treatment (N. K. Singh et al., 2023) and are increasingly being applied to treatment wetlands. These models can be trained on individual treatment wetland data (Yang et al., 2023) or data from wide ranges of treatment wetlands (Nguyen et al., 2022). Models trained on many different treatment wetlands enable the findings of individual studies to be combined in a model that uses common design parameters as the input variables (called features in machine learning). These models can be interpreted to expose their underlying logic (Molnar, 2024), potentially enabling insights such as those from simulations of the mechanistic models. 18 1.3 Problem statements Research gaps remain in studies of specific treatment wetland design, operation, and climatic instances, as well as how modelling is used for developing a broader understanding of their dynamics. 1.3.1 Limitations of temperature inclusion and model interpretation of treatment wetlands Treatment wetland machine learning models have not yet sufficiently addressed colder temperatures, nor provided an adequate explanation of how their models interpret the relationships within the data. Models trained on data from many publications have used air and water temperature interchangeably, which is an incorrect assumption in colder weather. To our knowledge, no machine learning models have been trained to predict wetland water temperature either, despite its importance in design. Finally, machine learning results are often limited to prediction accuracy, which does not progress the understanding of treatment wetlands, or are limited to showing the average relationship of a single algorithm, which does not show potential variability or uncertainty. 1.3.2 Batch operation of C. utriculata treatment wetlands through seasonal changes Studies of the wetland plant C. utriculata show better cold temperature performance than any other plant species but have been largely conducted with batch microcosms in unrealistic conditions, leaving a gap in understanding for more practical design. Long hydraulic retention times were used to emphasise the impact of the plants, made possible by continually adding water to compensate for evapotranspiration (W. C. Allen et al., 2002a). In real applications these retention times would dry out the treatment wetlands in summer, 19 potentially killing the plants or changing the dominant species. If retention time is shortened, the roots would have less time to deliver oxygen through their roots, potentially lowering their beneficial effect in more applicable conditions (Taylor et al., 2011). Neither evapotranspiration-caused water level reduction throughout each batch nor the impact of shorter batch times in colder temperatures have been studied to our knowledge. Finally, we do not know if cutting the above ground stalk would eliminate the cold-temperature benefit over unplanted beds. 1.4 Objectives General objective: develop an integrated understanding of the dynamic relationships between design, operational, and environmental factors affecting cold-climate treatment wetland performance through complementary machine learning analysis and controlled experimentation. Specific objectives: 1) To quantify relationships between temperature, design parameters, and treatment efficiency through interpretation of machine learning models trained on published treatment wetland data, with specific emphasis on factors influencing cold-climate performance. 2) To experimentally validate and extend machine learning insights by evaluating C. utriculata's treatment capabilities under more realistic batch operation conditions, including quantification of: a) Evapotranspiration effects on bed hydraulics and treatment efficiency b) Cold-temperature organic matter removal mechanisms 20 c) Post-harvest treatment resilience in cold conditions 1.5 Thesis organization This thesis is composed of four chapters. Chapter one presented background theory, problem statements, and the objectives of this thesis. Chapter two presents the results of collecting published data on cold climate constructed wetland dynamics and training and interpreting machine learning models which predict effluent temperature and quality. Chapter three discusses insights and research gaps generated from the machine learning research, and how they can be addressed using microcosm experiments. Chapter four presents my experiments with batch constructed wetland microcosms using the plant species C. utriculata and discusses the results. Chapter five concludes the thesis and discusses potential future research directions. 1.6 Scientific outcomes and outputs Presenting various stages of this work while it was in progress allowed sharing and receiving information with the wider scientific community. These presentations included: • UNBC - Three-minute Thesis 2023: Wetland designs for wastewater treatment in cold climates • UNBC - Three-minute Thesis 2024: Treating wastewater through seasonally cold climates • UNBC - Faculty of Science and Engineering Poster Competition 2024: Planting Sedge to Improve Wastewater Treatment at Cold Temperatures • International Water Association - World Water Congress 2024: Seasonal Performance Dynamics for Subsurface Treatment Wetlands in Cold Temperatures 21 I also intend to publish this masters work through two manuscripts: • Manuscript 1: Interpretable Machine Learning Reveals Cold Climate Dynamics of Treatment Wetlands • Manuscript 2: Carex utriculata Improves Treatment in Three-Day Batch Operated Treatment Wetlands in Warm, Cold, And Harvested Conditions 22 Chapter 2 Interpretable Machine Learning Reveals Cold Climate Dynamics of Treatment Wetlands 2.1 Introduction Subsurface flow treatment wetland adoption in cold climate regions (mean temperature <0°C during coldest month) (Cui et al., 2021) has been hindered by the potential for temperaturedriven freezing and performance reduction (M. Wang et al., 2017). There are successful strategies in cold climates (Grebenshchykova et al., 2020; Kadlec & Wallace, 2009), but additional research is required for robust designs which combine previously researched strategies (Ji et al., 2020). Comparing and combining the results from cold climate constructed wetlands can be complex due to differing temperature and design factors. In colder conditions water temperature is primarily determined by the influent temperature (Muñoz et al., 2006) which has been rarely reported or considered in constructed wetland studies (Grebenshchykova et al., 2023). Cold-climate studies reporting similar air temperatures may experience different water temperatures. Furthermore, the effect of water temperature on microbial activity is a function of both the microbial community and the availability of substrates (Nedwell, 1999), which vary by design. Both water temperature and characteristics of the wetland design and operation must be known and accounted for to compare study results. The impact of design and operational variables have been previously investigated through supervised machine learning algorithms (Lam et al., 2024; Nguyen et al., 2022). In contrast to more traditional modelling techniques such as the P-K-C* model developed by Kadlec and Wallace (2009) or more complex mechanistic models such as BIOPORE (Samsó & Garcia, 23 2013), machine learning models can be trained with variables available prior to construction, rather than those measured or calculated afterwards. Machine learning models of constructed wetlands have also shown greater accuracy than more basic linear regression methods (Nguyen et al., 2022, 2024a). Despite their benefits, there are issues with the transparency of machine learning models which need to be addressed using interpretable machine learning methods (Hassija et al., 2024). As machine learning models are often black boxes, the lack of transparency in their logic necessitates distrust of the results (Moraffah et al., 2020). Interpretable machine learning is an active field of development which allows the models understanding of the relationships and importance of variables to be visualized (Molnar, 2024). These are valuable as model relationships do not necessarily match causation (Molnar et al., 2022), and equally accurate models can form different relationships (Fisher et al., 2019). To increase the robustness and value of modelling results, multiple models can be compared to each other (Gunasekaran et al., 2024) and established theory. Significant gaps remain in the use of machine learning for understanding how temperature, design, and operational dynamics affect constructed wetlands. The only machine learning studies found which used water temperature as a predictive variable (as opposed to air or a mixture of the two) were limited to data from individual studies (Akratos & Tsihrintzis, 2007a; Nivala et al., 2019). To our knowledge machine learning has never been used for predicting constructed wetland effluent temperature, despite the importance of knowing water temperature for design (Nedwell, 1999). Finally, interpretation of machine learning models has so far been limited. Algorithms have only been compared using their prediction performance (Yang et al., 2023), and interpretation plots have been limited to single 24 algorithms (Lam et al., 2024; Nguyen et al., 2022). This work fills the existing gaps in literature by: (i) using water temperature data from across published literature as a predictive variable for treatment performance, (ii) training models to predict effluent temperature, and (iii) comparing the interpretation of multiple algorithms to each other and the literature. By filling these gaps, we take another step towards understanding how temperature, design, and operational dynamics affect the performance of treatment wetlands, and how machine learning can further this goal. 2.2 Methods 2.2.1 Data collection Data was obtained from peer reviewed published articles; both directly from the published papers as well as wider or more detailed datasets shared directly by the authors of those articles where available. To be included, data was required to be from subsurface treatment wetland systems and have information on features considered critical based on previous studies. Features are any individual variables, characteristics, or properties of a system used by a machine learning model to make predictions. Critical features were media diameter, media depth, plant status (unplanted, planted, harvested), flow direction (horizontal, vertical), saturation (unsaturated/free draining, permanently saturated, tidal/alternating/reciprocating, aerated, and single-fill-and-drain), loading interval (hours between intermittent doses), hydraulic loading rate (HLR, L/m2/d), and effluent temperature (Kadlec & Wallace, 2009; Langergraber et al., 2020). Where systems were small scale, not buried, and did not report excessively warm influent the effluent temperature was assumed to be equal to the air temperature due to the high rate of heat transfer (Ouellet-Plamondon et al., 2006). 25 Data was split into temperature and treatment datasets to train machine learning models to predict effluent temperature and treatment performance respectively. For the treatment dataset, measurements of organic matter and ammonia influent and effluent were also required. For the temperature dataset, systems had to be outside buried systems (OuelletPlamondon et al., 2006) and have data on air and influent temperature (Prost-Boucle et al., 2015a). In the few cases where air temperatures were not reported, these were obtained where possible from historical weather data for the location and time period of the study. In cases of multi-stage constructed wetland systems, each stage was considered a separate system, and data collection was further limited to stages where influent and effluent data were available for each stage. All temporal variables were required to be reported as individual data points or be the mean or median over a maximum of one season in order to preserve meaningful relationships between variables. Non-critical values shown to have theoretical importance in literature were also collected where available, but their absence did not exclude new data sources from being added. These included anodes (W. Liu et al., 2022), rest time between batch operations (Jia et al., 2010), and insulative properties such as depth and presence of snow (Grebenshchykova et al., 2020), detritus, and insulation (Kadlec & Wallace, 2009). Unbalanced, high cardinality, correlated, or dependant features were removed from the resultant datasets where possible. Highly imbalanced design modifications shown to have a large impact on the target variable (e.g. anodes) prompted exclusion of treatment wetlands with the minority category. Highly correlated features were identified using Spearman’s rank correlation matrix (Liebetrau, 1983) and minimised to prevent obscuring the interpretability of the model (Farrar & Glauber, 1967; Molnar et al., 2022). Correlated features were kept if 26 both features had unique and potentially theoretically valuable information. In these cases, potential importance misattribution issues were subsequently investigated using accumulated local effect plots after model generation. 2.2.2 Model construction Data was processed in R and modeled using the R tidymodels framework (Kuhn & Wickham, 2020). Text-based factors were converted to dummy variables, all values were normalized, and a near-zero variance filter was applied to ensure highly unbalanced or sparse features did not exist. Data was divided into training (80%) and test (20%) datasets as done with similar studies (Nguyen et al., 2024b), stratified by 5°C increments of effluent temperature to ensure equal representation of temperature ranges within the training and test datasets. Five-fold stratified cross-validation was used for tuning to reduce the risk of overfitting while maximizing the use of all available data (S. Singh et al., 2022a). Post-cross validation, the best tuning parameter combination was automatically selected based on the optimal RMSE value across a twenty-point grid of tuning values. Models used were chosen according to their success in other publications which modelled subsurface constructed wetland performance (Nguyen et al., 2024a; S. Singh et al., 2022b), with the addition of a gradient boosting model (XGBoost) not yet tested for treatment wetlands which has outperformed other models for modelling denitrification of sludge (S. Liu et al., 2025). Non-linear models were used exclusively due to the established non-linear nature of relationships between design variables and both treatment rate and effluent temperature in constructed wetlands (Kadlec & Wallace, 2009). Boosted trees via the extreme gradient boosting (XGBoost) R package (GBM) is a sequential tree building model for sparse data that excels at complex, nonlinear interactions and attempts to generalize and reduce 27 overfitting (Chen & Guestrin, 2016; Sarker, 2021). Random forests via the Ranger R package (RF) aggregates multiple decision trees that averages independent predictions, increasing robustness but reducing sensitivity to more subtle relationships (Wright & Ziegler, 2017; R. Xu et al., 2016). Polynomial support vector machines with the radial basis function via the kernlab R package (SVM) transform data to a higher dimensional space, where it forms a hyperplane that can capture complex non-linear relationships (Bonaccorso, 2018; Karatzoglou et al., 2004) but can perform less well for noisy data (Sarker, 2021). Cubist via the cubist_rules R package (Cubist) is a rule-based tree model followed by regression models to increase flexibility. Each model brings a different approach to the table, increasing the value of comparing each model. 2.2.3 Interpretation To determine the impact of each feature within each model and algorithm, SHAP (SHapley Additive exPlanations) values were generated using the fastshap package in R (Greenwell, 2024). To provide a more robust assessment of feature importance, the median SHAP value across each of the four algorithms was plotted (Fisher et al., 2019). Key strong interactions between features were identified using Friedman’s H-statistic method (Friedman & Popescu, 2008). The H-statistic method measures the average fraction of variance explained by interaction with a feature and was calculated using the iml R package (Molnar et al., 2018). Interactions with values above 0.1 and/or standard deviations below 0.1 were targeted for interpretation. To visualize the impact of features which were not inherently correlated with other features, partial dependence plots (PDPs) and individual conditional expectation (ICE) plots were used. ICE curves each show how different individual random combinations of feature inputs 28 would respond to varying a chosen feature, with each curve plotted as a thin partially transparent line. PDPs, shown with bold thicker lines, show the average marginal impact of varying a chosen feature across all feature combination possibilities (Friedman, 2001). To visualize the impact of features inherently correlated with other features, accumulated local effects (ALE) plots were used. ALE plots only consider feature combinations within specific intervals of the feature(s) being considered to reduce misleading results that could be introduced from unrealistic combinations of inherently correlated features. The effects in each local interval are then calculated as differences from the mean prediction, and connected to form a curve (Apley & Zhu, 2020). ALE plots show either the main effect or the interaction effect (Molnar, 2024), not the combined effect, therefore main and interactive effects were both plotted where applicable to show the full picture. ICE, PD, and ALE plots were used to compare each models understanding of relationships between features to the findings in literature. Studies which study the average effect of a variable across different design permutations are limited. The P-K-C* model (Kadlec & Wallace, 2009) is one, which was compared to our HLR and influent strength findings since these are both variables within the P-K-C* equation. Remaining P-K-C* variables were chosen based on mean HLR and organic matter influent strength in our datasets, along with mean C*, P, and k values as reported by Kadlec and Wallace for our mean influent strength. 29 Table 2-1: P-K-C* variables which correspond to the mean treatment wetland features in our machine learning data sets. These are used to provide predictions of effluent concentrations to compare with model interpretation plots. Variable C* (background concentration) P (tanks in series constant) k [m/y] (areal removal rate constant) HLR [L/m2/d] Vertical unsaturated 2 Horizontal saturated 10 2 3 187 25 90 27 2.3 Temperature results and discussion 2.3.1 Data collection Data collection criteria are summarised in Table 2-2, resulting in 7 publications, 27 unique treatment wetlands, and 1933 data points (see list of papers in Table A-1). Influent temperature was the least reported critical feature, significantly limiting data which datasets could be included. Radiation, ground transfer, and evapotranspiration components of the energy balance equation (Kadlec & Wallace, 2009) were not directly included as critical features due to lack of available data. These may reduce accuracy of summer temperatures but are considered negligible for prediction of temperatures in the winter due to snow and ice cover (Grebenshchykova et al., 2023). Insulating features suffered from high cardinality and imbalance. Insulative materials in the collected data included mulch, plastic, foam, glass beads, and gravel while collectively representing only 2% of datapoints. Presence and depth of snow was not consistently reported, especially as winter air temperature variation can cause periodic melting (Grebenshchykova et al., 2020; Muñoz et al., 2006) that requires frequent reporting of snow presence and depth. Plant detritus thickness was variable but not consistently reported. To lower imbalance and cardinality while still accounting for insulative materials, a single 30 numerical R value [m2K/W] feature was calculated for each system based on reported insulation thickness and estimated thermal properties (Kadlec & Wallace, 2009). Table 2-2: Features considered for training the temperature model and the reasons for including or not including them. Feature Influent temperature Air temperature Saturation HLR Insulation Plant status Snow Ground temperature Evapotranspiration Radiation Imbalanced variable 1 Expected Importance High (Grebenshchykova et al., 2023) Moderate (Prost-Boucle et al., 2015a) Moderate (Prost-Boucle et al., 2015a) Moderate (Grebenshchykova et al., 2023) Moderate (Kadlec & Wallace, 2009) Moderate (Agnieszka & Anna, 2013) Moderate (Grebenshchykova et al., 2023) Low (L. Zhang et al., 2015) Low2 Low2 Availability Low High High High Moderate1 Moderate Low Very low Low Very low Included Yes Yes Yes Yes Yes No No No No No Low only in cold-temperature conditions, which are the focus of this study. High in warmer 2 temperatures (Wallace & Nivala, 2005). 2.3.2 Model performance XGBoost achieved the highest accuracy, however similar performance to Cubist and RF shows that comparing the interpretation results of each model is necessary to understand their differences (Figure 2-1). The good accuracy from XGBoost aligns with its strategy of continual sequential tree improvement. Despite lower fit metrics, SVM was retained in the analysis as its less flexible nature provides a valuable comparison during model interpretation exercises (Awad & Khanna, 2015). To the best of our knowledge there are no other studies which use machine learning techniques to predict temperatures in constructed wetlands. Empirical models exist with R2 values of 0.86 to 0.94 for horizontal subsurface constructed wetlands (Kadlec & Wallace, 2009) but rely on inputs of annual average water temperature and amplitude of the sinusoidal 31 annual pattern. Our model has the advantage of using values obtainable prior to construction as well as achieving higher fit metrics. Figure 2-1: Effluent temperature predictions plotted versus actual effluent temperature for each model algorithm tested, as well as their R2 and RSME goodness of fit scores. 2.3.3 Model interpretation The SHAP values for each feature were similar between models, however model disagreement increased for interaction strengths and increased again for specific interactions (Figure 2-2). The SVM model was the largest outlier in both interaction plots, indicating differences in modelled relationships compared to the other models. Omitting SVM results, interaction strengths shows that interactions were responsible for no more than 20% of the variance in predicted effluent temperature per feature, and main effects dominate predicted effluent temperature. Interaction strengths of influent temperature, delta, and HLR were highest, and the strength of specific interactions between these features which also had lower standard deviations were delta:HLR and delta:influent temperature. 32 Figure 2-2: Ranked importance of features and interactions between features for the prediction of effluent temperature. SHAP values show the overall importance of each feature. Interaction values show how much variance in predictions is attributed to interactions of either a particular feature (interaction strength) or a specific interaction (shown with “:” between the interacting features). The SHAP values in Figure 2-2 agree with a study of vertical flow filters (Grebenshchykova et al., 2023) which showed that influent temperature has the highest influence. The same study concluded that HLR was more important than air temperature, though air temperature was noted to be especially important at the beginning of the cold season. Our use of delta inherently accounts for the larger difference between air and influent temperature during the seasonal temperature shift, which may account for delta having higher importance than HLR. Further investigations of feature relationships using ALE plots (Figure 2-3) show linear effects of influent temperature and delta between influent temperatures of 6 and 20°C and -15 and 2. These linear relationships align with those in the temperature balance equation (R. Kadlec & Wallace, 2009). The ALE plots show less linear relationships in higher and lower delta and influent temperature conditions. In warm conditions this may be accounted for by evapotranspiration consumption of energy, which is the largest warm temperature energy loss 33 (R. H. Kadlec, 2006) that was not reported enough to be included. In colder conditions there is ice and snow formation, which slows the effect of delta (Wallace & Nivala, 2005) and also did not have sufficient data to be included. Figure 2-3: Accumulated local effect of each feature in the model. The marks along the x-axis indicate data point presence, giving an approximation of data distribution. ALE plots of the interactive effect between HLR and delta (Figure 2-4) showed that HLRs less than 10 – 15 L/m2/d cause delta to have a larger (0 – 1°C) impact on predicted effluent temperature at both higher (10) and lower (-30) delta values. This is supported by the temperature balance equation (Equation 1), with the same effect observed by other studies (Grebenshchykova et al., 2023) at higher HLRs (70 – 140 L/m2/d) than our ALE plots had clear patterns for. This suggests that there is a wider effect than was clearly captured in our model interpretation plots. 34 Figure 2-4: The accumulated local effect of delta on the predicted effluent temperature depending on the specific algorithm (right), and the accumulated local interaction effect between delta and HLR for the Cubist (top-left) and XGBoost (bottom-left) algorithms. The colour legend for the interaction effects is shown in the top-right. Presence of training data points are shown as black dots as well as black lines along the x and y axis. No effect of saturation was observed in the original model. As saturation type has been previously shown to effect temperatures within the bed (Prost-Boucle et al., 2015a), a simplified model was trained without the HLR, as it is highly correlated with saturation type, to see how model interpretation changed. The simplified models disagreed on the magnitude of the saturation effect, but each algorithm showed alternating systems to be most affected by delta, followed by aerated, percolating, and lastly saturated systems (Figure 2-5). This matches the order of both the rate of air inflow and the mean HLR from largest to smallest. As increased airflow should increase the influence of delta and increased HLR should lower the influence of delta, this suggests that saturation differences may have a larger impact than HLR and deserves further study. 35 Figure 2-5: ALE of delta on the predicted effluent temperature for different algorithms (separate plots) and saturation types (different coloured lines). Training data presence is noted by marks along the x and y axis, coloured by saturation type. ALE interactions between plants or insulation and delta or influent temperature showed only very slight (<0.25°C) effects on the predicted effluent temperature. Though marginal, effects showed planted and insulated systems warmer at cold temperatures, and colder at warm temperatures. Plants have shown to have a similar impact on bed temperature of unsaturated systems as sludge buildup in French vertical flow systems, greatly reducing the decline between doses even when harvested and removed prior to winter (Torrens et al., 2009). As this lowered bed temperature does not directly correspond with reductions in the effluent temperature (Prost-Boucle et al., 2015a) predicted by our model, this may account for the lowered effect. All specific interactions (Figure 2-2) with strengths above 0.1 or standard deviations below 0.1 were assessed, though not all are displayed here. Most interactions were both small and slightly different in each model, resulting in no consistent relationship. Interactions with large standard deviations often predicted large effects in values in regions of lower data density, and low effects within higher data density areas, indicating a lack of any robust relationship. 36 2.4 Treatment results and discussion 2.4.1 Data collection Data collection criteria are shown in Table 2-3, resulting in an initial 24 publications, 108 unique treatment wetlands, and 2350 data points (see detailed list in Table A-2). Effluent temperature was reported rarely in literature (with air temperature more common), and the frequency of feature availability often did not meet requirements, making these two restrictions the primary limiting parameters for data inclusion. Features recorded during data collection but with insufficient data in the resulting dataset to be included were organic nitrogen (available in 55% of total data), rest times between fill-and-drain batches (available in 70% of fill-and-drain data), and evapotranspiration (reported in 43% of total data points). High correlations eliminated dosing method and flow direction, which had high correlations with each other and saturation. All but one fill-and-drain system was batch dosed, all but two saturated systems were continuously dosed, and all vertical unsaturated systems were intermittently dosed. For flow direction, fill-and-drain systems were the only saturation type which had substantial numbers of systems in both vertical and horizontal configurations. Imbalanced data included recirculation and plant species. Recirculating systems represented only 2% of the final datapoints but were retained due to being identified as a solution for cold temperature systems (Ji et al., 2020) and having clear effects in the datasets which were available (Arias et al., 2022). Plant species was omitted due to imbalance and high cardinality, with Phragmites australis representing 80% of the data, with the remaining 20% composed of 23 other species. 37 Large high-resolution datasets made up 70% of the data points for treatment wetlands older than 1 year, resulting in treatment wetland age serving as a proxy for seasonal temperature changes. To account for age while preventing its use as a proxy for season, ages over 1 year were rounded to one significant figure. Table 2-3: Features considered for inclusion in the treatment model training data, including whether they were included and the rationale. Feature COD sCOD BOD5 CBOD5 NH4-N Org-N Influent temperature Bed temperature Effluent temperature Air temperature Saturation Loading interval Batch rest time Dosing method Flow direction Hydraulic loading rate Recirculation rate Media diameter (smallest) Media diameter (median to largest) Influent distribution system Plant status Plant species Plant family Evapotranspiration Age Expected Importance Target variable Target variable Target variable Target variable Target variable Moderate (Thorén et al., 2004b) High (Pucher & Langergraber, 2019) Low (Prost-Boucle et al., 2015b) High (Pucher & Langergraber, 2019) Low (Prost-Boucle et al., 2015b) High (Nivala et al., 2013) Moderate (Langergraber et al., 2020) Moderate (Jia et al., 2010) Uncertain, high correlation with saturation Uncertain, high correlation with saturation High (Pucher & Langergraber, 2019) Availability High Low Very low Moderate Moderate Low Very low Very low Low High High Moderate Low High Included Yes Yes Yes Yes Yes No No No Yes No Yes Yes No No High No High Yes Moderate (Foladori et al., 2013) High (Pucher & Langergraber, 2019) Low High Yes Yes Low (Langergraber et al., 2020) High No Moderate (Pucher & Langergraber, 2019) Low (Taylor et al., 2011) High (Taylor et al., 2011) Low No High Imbalanced, high cardinality Imbalanced Low Yes No High Yes Moderate (Taylor et al., 2011) Moderate (Ouellet-Plamondon et al., 2006) Moderate (Ding et al., 2021) No No 38 2.4.2 Preliminary modelling and data division Preliminary modelling for prediction of organic matter used all data in the treatment dataset but resulted in misattribution issues. Feature importance scores (not shown) gave depth the highest median importance across the models, while saturation had low importance. Similarly, ALE plots showed that the models predicted lower effluent concentration depending on depth and HLR rather than the saturation type. This is contrary to the established understanding, which designates the type of wetland as the largest factor for treatment due to oxygen transfer differences (Nivala et al., 2013). Correlation matrices showed that depth, HLR, saturation, organic matter, and ammonia features were highly correlated due to the typical performance and design parameters of each type of wetlands (Table 2-4), likely causing misattribution issues in the models. Removing HLR and depth from the preliminary model resulted in saturation having the highest SHAP value, with ALE plots showing large interactions between saturation and other features. This confirmed collinearity issues and showed that the effect of other features was highly dependant on saturation. Similar machine learning studies which have used vertical and horizontal flow data in the same model have also found HLR and filter depth to have high importance (Lam et al., 2024; Nguyen et al., 2022). These correlation issues may have been the cause, with Lam showing PDPs for HLR and depth with a similar step pattern to those generated by all-inclusive models of this dataset. The lower step corresponds with the typical design parameters of saturated systems, while the higher step corresponds with those of unsaturated systems. Based on our results, future studies should avoid correlation issues by separating saturation 39 types (Soti et al., 2023) or removing correlated design features and retaining only system type (S. Liu et al., 2025). Table 2-4: Mean values of features correlated with saturation type in the treatment data set. as well as how many unique treatment wetlands were collected within each saturation type. Saturation Systems HLR [L/m2/d] Media Depth [m] aerated 9 108 ± 20 0.9 ± 0.1 alternating 1 169 ± 13 0.9 ± 0 fill-and-drain 39 18 ± 25 0.3 ± 0.2 percolating 22 128 ± 117 0.7 ± 0.2 saturated 34 30 ± 28 0.4 ± 0.1 Influent Effluent COD: 432 ± 140 CBOD: 298 ± 83 NH3: 65 ± 17 CBOD: 277 ± 53 NH3: 55 ± 15 COD: 369 ± 177 BOD5: 173 ± 50 NH3: 17 ± 11 COD: 373 ± 215 CBOD: 269 ± 65 BOD5: 131 ± 134 NH3: 43 ± 22 COD: 513 ± 194 CBOD: 294 ± 89 BOD5: 98 ± 49 NH3: 57 ± 26 COD: 24 ± 8 CBOD: 2 ± 3 NH3: 1 ± 1 CBOD: 4 ± 4 NH3: 8 ± 10 COD: 61 ± 51 BOD5: 21 ± 16 NH3: 17 ± 7 COD: 24 ± 10 CBOD: 12 ± 17 BOD5: 18 ± 32 NH3: 8 ± 11 COD: 63 ± 51 CBOD: 54 ± 23 BOD5: 30 ± 22 NH3: 45 ± 26 Subsequent modelling used separated datasets, which enabled retention of HLR and depth while increasing interpretability. The dataset was divided by saturation to eliminate high correlations, with batch and saturated systems grouped due to similar design parameters and performance (Table 2-4). As saturation was shown to heavily influence the effect of other features, dividing the dataset also served to account for saturation in each interpretation method. This allowed for more useful interpretation of interactions since there are no interpretation techniques for interactions of more than two features (Molnar, 2024). As a result of the data division, aerated and alternating systems were removed from modelling due to high treatment consistency and data limitations (Table 2-4). These systems had high and consistent treatment rates across all designs and conditions reported, reducing 40 the benefit of modelling for the data distributions available. In addition, only one alternating system had data which met requirements and 99% of the aerated datapoints were from the same paper, providing less opportunity for insights not already captured by the original publications. The decision was therefore made to remove them from the modelling analysis, and the final data used for modelling is shown in Table 2-5. This final data totalled 24 publications, 100 unique treatment wetlands, and 1693 data points. Table 2-5: The range, mean, and standard deviations of numeric features in the separated treatment datasets. Organic matter Saturated Unsaturated min – max vertical Feature Organic matter effluent [mg/L] Organic matter influent [mg/L] NH3 effluent [mg/L] NH3 influent [mg/L] Minimum media diameter [mm] HLR [L/m2/d] Depth [m] Age [y] Effluent temperature [°C] Recirculation rate [%] Loading interval [h] (mean±sd) min – max (mean±sd) 0 - 300 (61±39) 59 - 900 (370±160) 0 - 140 (16±20) 110 - 550 (300±100) N/A N/A N/A N/A 0.25 - 30 (7.7±6) 3.3 - 380 (27±25) 0.25 - 1.1 (0.42±0.11) 0.05 - 9 (2.6±2) 0.1 - 27 (12±5.7) 0 – 300% (0.62±11) 0.06 - 4 (1.7±1.4) 22 – 250 (90±28) 0.1 - 1 (0.78±0.18) 0.1 - 6 (1.3±0.91) 3.1 - 27 (13±5.4) 0 – 200% (0.023±0.16) 1 - 12 (2.1±1.9) N/A NH3 Saturated Unsaturated min – max vertical (mean±sd) min – max (mean±sd) N/A N/A 0 - 900 (350±180) 0 - 160 (44±27) 0.57 - 270 (54±24) 0.25 - 30 (8.2±6) 5.5 - 380 (27±26) 0.25 - 0.7 (0.43±0.12) 0.05 - 6 (2.3±1.6) 0.1 - 27 (13±5.7) 110 – 550 (290±100) 0 - 54 (11±12) 14 - 150 (54±19) 0.06 - 4 (1.7±1.4) 22 – 250 (90±29) 0.1 – 1 (0.78±0.18) 0.1 - 6 (1.3±0.92) 3.1 - 27 (13±5.5) 0 – 200% (2.3±16) 1 - 12 (2.1±2) N/A N/A 41 2.4.3 Model performance Model fit primarily differed by dataset rather than algorithm (Table 2-6). In general, our models demonstrated similar or better accuracy than those reported in previous machine learning studies of treatment wetlands. However, the prediction of organic matter in saturated systems showed a weaker fit. For instance, horizontal saturated treatment wetlands modeled with SVM have achieved an R2 of 0.4 (S. Singh et al., 2022b). In another study, multiple treatment wetland types were modeled using Cubist, RF, and SVM, finding that Cubist performed best with an R2 of 0.84 for both ammonia and organic matter models (Nguyen et al., 2022). Our fit metrics (Table 2-6) were lower for organic matter but corresponded to the standard deviations of our dataset (Table 2-5). The relatively poorer prediction for saturated systems may be attributed to the wider range of effluent concentrations in our dataset (61±39 mg/L organic matter), compared to Singh’s study which reported (10±12 mg/L BOD, 42±14 mg/L COD). Table 2-6: Fit metrics for each model and algorithm that is predicting effluent concentrations. Dataset XGBoost rmse / R2 Cubist rmse / R2 SVM rmse / R2 RF rmse / R2 saturated 37 / 0.33 38 / 0.31 44 / 0.21 37 / 0.31 vertical unsaturated 16 / 0.58 11 / 0.81 21 / 0.28 13 / 0.71 Ammonia saturated 6.8 / 0.93 8.4 / 0.90 9.3 / 0.87 7.7 / 0.91 Ammonia vertical unsaturated 2.9 / 0.92 4.0 / 0.87 6.5 / 0.68 3.7 / 0.89 Parameter Organic matter Organic matter 2.4.4 Model interpretation SHAP values were highly dependant on both the saturation and whether ammonia or organic matter was being predicted (Figure 2-6). Feature importance scores have previously been 42 shown for machine learning of treatment wetlands (S.-Z. Zhang & Jiang, 2024), but with different features, and comparison of machine learning feature importance scores across saturation types and algorithms has not been done to our knowledge. In our model the large variability between models shows that feature importance scores of a single model using data from all saturation types would have limited value. Comparison of algorithms shows larger standard deviations for features with higher median SHAP values, but similar feature rankings, giving higher confidence to the rankings. Interaction importance scores were not generated for the treatment models due to the processing power required and visualization difficulty for an increased number of models and features. Figure 2-6: SHAP values show the relative importance of each feature depending on if ammonia or organic matter is being predicted (top labels, first line) and if the training data is from saturated or unsaturated treatment wetlands (top labels, second line). Coloured dots show SHAP values of specific algorithms, while the box plots show their median importance and distribution. 2.4.4.1 Influent loading Organic matter effluent prediction showed a lower response to organic matter influent concentrations than modelled using P-K-C*. Unsaturated models predicted an average 2 mg/L rise in effluent per 100 mg/L of influent increase, while saturated predicted an average 5 mg/L rise. As comparison, the P-K-C* equation (Kadlec & Wallace, 2009) modeled rises of 20 mg/L for vertical flow and 40 mg/L for horizontal flow. This is representative of the high removal rates in our data. As our model was trained on individual points rather than median 43 values, it may also be a consequence of robustness (Seintos et al., 2024; Sultana et al., 2016) to temporary changes in influent concentration. Figure 2-7: PDP and ICE plots show how organic matter influent (x-axis) [mg/L] changes the predicted organic matter or ammonia effluent (y-axis) depending on the algorithm (colored lines) and saturation type (separate plots). The step response in the unsaturated model’s predicted ammonia effluent as a function of influent strength (Figure 2-7) was consistent between algorithms. The response also aligns with slow nitrification reported prior to BOD being reduced below 20 mg/L, due to competition for oxygen between nitrifiers and heterotrophs (Headley et al., 2005). For saturated treatment wetland models, low removal of ammonia was predicted regardless of the influent organic matter concentration, indicating that it is not the limiting factor. This low removal rate is shown more clearly in Figure 2-8, where influent and effluent concentrations of ammonia are highly correlated for both saturated and unsaturated treatment wetlands. 44 Figure 2-8: PDP and ICE lines show how ammonia influent [mg/L] influences the predicted ammonia effluent for saturated (left) and unsaturated (right) treatment wetlands depending on the algorithm (coloured line). Unsaturated vertical wetlands had no training data with influent ammonia greater than 100 mg/L. The effect of HLR in our models (Figure 2-9) was greatly diminished compared to literature. Between 0 and 100 L/m2/d the P-K-C* model showed organic matter effluent rising 400 mg/L for saturated systems and 200 mg/L for unsaturated systems. In our models only the saturated model predicted organic matter increased, and it was an order of magnitude lower. It is possible that the lowered effect of HLR in the models was due to strong buffering capacities (Weerakoon et al., 2016), as we used multiple data points from each system, some of which varied in HLR. 45 Figure 2-9: PDP and ICE plots of how HLR [L/m2/d] influences each algorithms predicted organic matter or ammonia effluent, depending on the saturation of the treatment wetland. 2.4.4.2 Physical design parameters Media diameter interpretation plots generally showed that smaller media can achieve lower effluent concentrations for both saturated and unsaturated systems. This is well understood for unsaturated vertical systems, for which simulations calibrated on experimental data have been published (Langergraber et al., 2020). For both ammonia-N and organic matter effluent, Langergraber’s results show the largest increase between 0.06 mm and 1 mm, with decreased benefit above 1 mm. While our overall organic matter increase from 0 to 4 mm is similar, our decreases between 0.06 mm and 1 mm and subsequent jumps in concentration are not aligned with the theoretical understanding. Ammonia-N results were more consistent between algorithms but did not show the larger difference in effluent concentrations between 0 and 1 mm that simulations predict. 46 Figure 2-10: PDP and ICE plots showing how the minimum media diameter [mm] influences the predicted ammonia and organic matter concentrations for each algorithm, depending on the saturation of the treatment wetland. For our saturated models, media size had a 10 – 30 mg/L impact, despite uncertainty in publications. While some studies showed smaller media slightly increasing organic matter removal (He & Mankin, 2007), others have not (Akratos & Tsihrintzis, 2007a). Likewise, while total Kjeldahl nitrogen (TKN) has shown to be reduced by smaller media (Akratos & Tsihrintzis, 2007a), the same conclusions do not appear to be present for ammonia-N. Our models show a clear effect on organic matter concentration, and a small effect for ammoniaN. The increased effect at smaller media sizes that other studies showed for unsaturated systems was predicted here for saturated. If this is the real relationship, it may have been better distinguished here due to the higher number of systems (88) in the saturated model compared to the unsaturated (22). The effect of depth in our saturated models (Figure 2-11) aligned with the clear impact in literature, with shallower systems performing better (Aguirre et al., 2005). For organic matter, changes below 0.6 m had increased effect for all but the SVM model. This would 47 correspond with the quicker percentage increase of bed depth near the water surface where oxygen transfer occurs. Unsaturated models showed the reverse effect, with removal increasing with depth. Unsaturated systems allow air to flow into their full depth between loadings, preventing increased depth from causing anaerobic conditions. Increased filter depth has shown to benefit all treatment parameters (Torrens et al., 2009), however this may be contingent on installation of passive aeration tubes which allow airflow to the bottom of the bed, better enabling aerobic conditions (Stefanakis & Tsihrintzis, 2012). Presence of aeration tubes was not included in the model to reduce dimensionality but could be added in future modelling efforts. Figure 2-11: PDP and ICE plots showing how the depth of the treatment wetland bed [m] influences the predicted ammonia and organic matter concentrations for each algorithm, depending on the saturation of the treatment wetland. 2.4.4.3 Temperature Responses to temperature varied significantly by saturation and parameter (Figure 2-12). Model responses generally aligned with established theory; removal is more sensitive to 48 temperature in less aerobic systems (Nedwell, 1999), and nitrifying bacteria are more sensitive to temperature than heterotrophic bacteria (Salvetti et al., 2006; Varma et al., 2021). Despite their sensitivity, saturated system PDP curves showed no ammonia-N response to temperature due to lack of predicted removal even at warmer temperatures. ICE curves revealed that feature combinations with zero effluent ammonia-N at warm temperatures show increased concentrations as temperatures lower. ICE curves also showed alternate responses for the other plots, indicating interactions with other features. Figure 2-12: PDP and ICE plots showing how effluent temperature [°C] influences the predicted organic matter and ammonia effluent concentrations for each algorithm, depending on the saturation of the treatment wetland. Plant presence effected effluent differently dependant on both temperature and saturation of the treatment wetlands (Figure 2-13). For the general saturated model, planted systems were predicted to have significantly reduced sensitivity to temperature decline, aligning well with 20-day batch studies included in the training dataset (W. C. Allen et al., 2002b; Taylor et al., 2011). However, these results must be taken with caution, as the 20-day batch time (3 – 7 L/m2/d) of these 103 data points represents a very low loading rate, and many of the plant species had uniquely high cold temperature performance compared to the more common 49 Phragmites australis. In addition, a study of batch (20-day) and continuous flow (5-day) microcosms showed minimal or no effect of plants in cold continuous flow microcosms, in contrast to large effect in batch microcosms (Stein et al., 2003). While the discrepancy between contact times in Stein’s study is obviously large, this deserved further investigation using modelling. Figure 2-13: PDP plots showing how effluent temperature [°C] influences the predicted effluent organic matter concentrations for all saturated systems (continuous and batch), depending on the algorithm (from left to right), and plant status (coloured lines). A supplementary model was created with only continuous flow saturated treatment wetlands, omitting all batch-operated microcosms. This model resulted in Figure 2-14, which predicts a more subdued but persistent benefit of plants at colder temperatures. This continuous flow saturated data set has a higher HLR (30 ± 28 L/m2/d, Table 2-5) than that of the batch microcosms (3 – 7 L/m2/d). It also consists of 80% Phragmites australis, which has been shown to be less effective at improving treatment in colder temperatures than other species (Taylor et al., 2011). While the discrepancy between Figure 2-13 and Figure 2-14 may be partially due to the difference in plant species, it may also be due to plants being less effective at higher HLRs or with continuous flow operation. Figure 2-14 is a more realistic relationship for the bulk of treatment wetlands, as continuous flow wetlands planted with Phragmites australis are more common than batch or other species outside of microcosm studies, and the average HLR is more aligned with typical design recommendations (Dotro et 50 al., 2017). For areas where P. australis is invasive or batch wetland operation is of interest, more research is necessary. Figure 2-14: PDP plots showing how effluent temperature [°C] influences the predicted effluent organic matter concentrations for continuous flow saturated treatment wetlands (no batch), depending on the algorithm (from left to right), and plant status (coloured lines). For unsaturated models, only a small decrease in effluent concentrations was attributed to plants in cold conditions (Figure 2-15), and only for the Cubist algorithm. To our knowledge this effect has not been reported for unsaturated systems, giving more uncertainty to this effect. A lack of measurable effect would align with the literature. Unsaturated treatment wetlands obtain more oxygen from their operation than theoretically obtainable from plants, therefore plant supplied oxygen would be a smaller percentage of the total oxygen available to them (Nivala et al., 2013), and make a proportionally smaller difference than in saturated operation wetlands. In warm conditions unsaturated models predicted that plants cause higher organic matter concentrations (Figure 2-15). This reveals one of the causes of the unsaturated model’s unintuitive improvement of organic matter removal as temperatures decrease. Literature on the effect of plants in unsaturated systems is scarce. A slight overall improvement in organic matter removal was reported for planted unsaturated systems (Torrens et al., 2009), but was not distinguished by temperature range, preventing adequate comparison here. Effluent 51 concentration due to evapotranspiration (Nivala et al., 2019) is a possible explanation for the model prediction, but has not been observed to increase effluent concentrations for unsaturated systems from our review of literature. Figure 2-15: PDP plots showing how effluent temperature [°C] influences the predicted effluent organic matter concentrations for vertical unsaturated systems, depending on the algorithm (from left to right), and plant status (coloured lines). Models predicted consistently lower ammonia effluent for planted systems regardless of saturation or temperature (Figure 2-16). Plant uptake of ammonia occurs during plant growth (Vymazal, 2007) in warm conditions, while root oxygen release increases in colder conditions (W. C. Allen et al., 2002b), providing oxygen for denitrification. Figure 2-16 suggests that plant uptake has a larger impact than root oxygen-caused denitrification for saturated systems. For unsaturated the additional root oxygen appears to be more beneficial. Studies of unsaturated and saturated system ammonia removal at both warm and cold temperatures were not found for comparison. 52 Figure 2-16: PDP and ICE plots showing how effluent temperature [°C] influences the predicted effluent ammonia-N concentration, depending on the saturation (upper and lower plots), algorithm (from left to right), and plant status (coloured lines). A limitation to all plant related prediction plots is that the evapotranspiration rate (and therefore also root oxygen supply and ammonia uptake) is a function of the air temperature, solar radiation, growth stage, and plant species (Kadlec & Wallace, 2009), none of which are included in the model. The strongest experimentally demonstrated impacts of plants on organic matter removal (which were included in model training data) were shown to occur at a water and air temperature of 4°C, when the air temperature was cold enough to stop or slow plant evapotranspiration, but not cold enough to wither the plants (Taylor et al., 2011). Our collected data from published studies of outdoor buried systems shows that a water temperature of 4°C corresponds to air temperatures ranging from -15 to 11°C (Figure 2-17). There is a lower air temperature limit below which oxygen supply through plant stalks lowers or stops; P. australis for example was found to have insignificant root radial oxygen loss below air temperatures of 2.4°C (Q. Wang et al., 2015a). For ammonia, evapotranspiration 53 driven uptake will also slow and stop as air temperatures approach zero; the common Typha spp. for example, starts to turn yellow at 5°C (Yan & Xu, 2014). In addition, studies have shown that plant effects are highly dependant on the species, which was omitted due to imbalance and high cardinality issues. Due to these unmodelled effects of air temperature and species, the ammonia and organic matter removal of planted systems is expected to vary from the behaviour shown by the PDPs dependant on the air temperature. Figure 2-17: Plot of influent (yellow) and effluent (blue) water temperatures versus air temperature for all outdoor buried systems in the temperature dataset. 2.4.4.4 Age For saturated models, PDP and ICE showed a quick decrease of organic matter effluent in the first four months to a year, aligning with reported startup period durations (Ding et al., 2021; Zidan et al., 2015). The subsequent 15 – 25 mg/L decrease in predicted organic matter effluent has also been observed for the treatment of dairy (Tunçsiper et al., 2015) wastewater. The dataset only included seven saturated systems older than three years, therefore additional studies may increase the certainty of this relationship. 54 For ammonia and unsaturated vertical models there was no clear relationship, with either no effect or disagreement between models. Start up effects have been reported for unsaturated vertical flow systems, with BOD5, COD, and NH4 all more erratic and higher during the first year (Stefanakis & Tsihrintzis, 2012). For French vertical flow systems, long term improvement due to age has been demonstrated to be similar to what our saturated model shows (Paing et al., 2015). The lack of effect for the unsaturated models may be a consequence of low data density, with only three vertical unsaturated systems older than two years. It is also possible that non-French unsaturated systems do not improve over time, since Paing attributed it to sludge buildup that improved distribution. Instead improved distribution can be designed for (Pucher & Langergraber, 2019). Figure 2-18: PDP and ICE plots showing how the treatment wetland age [years] affects the predicted organic matter and ammonia-N effluent concentrations, depending on the saturation and algorithm. 2.5 Limitations Models and their interpretations are a result of the distribution and range of data used. Lessrepresented treatment wetland types may not be accurately portrayed or insufficiently represented in the global interpretation methods which are biased towards treatment wetlands 55 with more data points. Models also predict based on correlations, while we are most interested in causation (Molnar et al., 2022). Patterns developed by models may match causal relationships, but this cannot be assumed blindly. For this reason relationships within models must be compared to specific publications and established theory where possible to assess the likelihood of causality and application to theoretical understanding. The majority of the datasets used are from small scale systems, which are not always representative of pilot or full scale systems, especially for planted systems (Brisson et al., 2024). Direct extrapolation of results between systems of different scale is assumed for this work but should be investigated further in the future. Seasonal comparison of a cold-climate hybrid treatment wetland and mesocosms using C. utriculata showed similar effluent quality results but unique microbial communities, showing that care must be taken when comparing (Ayotte et al., 2024). As data availability from pilot or full-scale treatment wetlands increases, wetland size could be included as a model feature or a dedicated model for full scale treatment wetland could be created. Temperature and water quality are both affected by past conditions which are not captured in these models, as they have been trained on point-in-time data points. Substantial heat capacity dampens temperature changes (Kadlec & Wallace, 2009), and microbial communities take time to stabilize from environmental changes. For example, lags in growth observed in spring temperature shifts (Kim et al., 2022) may indicate that more sensitive bacteria are not able to compete and regain population share until temperature is above a certain value, leading to possible differences in performance between spring and fall despite equal temperatures. This applies to plants as well, which may be fully matured (in fall) or not yet grown (in spring) with equal effluent temperatures. Research on this problem has shown 56 that long short-term memory (LSTM) neural networks can be effective, but should be combined with smoothed data, and have only been implemented for single treatment wetlands thus far (Yang et al., 2023). ALE and PD plots cannot visualize interactions of three or more features, despite these relationships existing in the models. While these higher order interactions are often less impactful than main effects or two-way interactions, it is possible for significant effects to be missed because of this. 2.6 Conclusion Machine learning models were able to predict effluent temperature with a RMSE of 0.7 °C, organic matter with a RSME of 20 – 40 mg/L, and ammonia with a RSME of 5 – 14 mg/L. Saturated systems had considerably more variable effluent quality than unsaturated vertical systems, resulting in lower RSME scores while achieving similar R2 values. Interpretation of the machine learning models successfully allowed comparison to the established understanding and identification of issues in model logic. Visualizing the range of response possibilities using ICE curves added considerable value. Areas of misalignment between model algorithms signalled more uncertainty, but did not always correlate with divergences of the curves from established theory. No single algorithm was clearly most aligned with the established understanding. While the SVM model had lower fit metrics and was not able to capture the less linear relationships, it also resisted the overfitting tendencies of the other more flexible models. High correlations between saturation type, HLR, depth, and effluent concentration caused all algorithms to make incorrect conclusions in preliminary 57 treatment modelling. Models must be interpreted with an understanding of data distribution, correlations, and wetland theory. Interactive relationships showed increased differences between model algorithms as well as inconsistencies with established understanding. This is unsurprising due to the increased degrees of freedom but does limit the building of new knowledge since it is these more complex relationships that are less understood. Interaction strength scores of individual features more successfully aligned with expected importance. Machine learning driven modelling can successfully extrapolate across results from multiple studies while accounting for differences in system designs. By leveraging model interpretation methods in combination with an understanding of the underlying theory of wetland systems, the average and range of effects of individual variables were effectively visualized. In this case, interpretation tools have provided a valuable review of the effects of not only temperature, but other key parameters. Using multiple modelling techniques in this analysis helped indicate where relationships were less certain and multiple explanations are possible, and where the relationship is more defined. The results of these methods can effectively communicate variable effects, allowing new insights into future system design and operation. These results must be interpreted with caution however, as modelling is based on correlation, not theoretical causation, and misattribution of cause can lead to incorrect predictions or insights. 58 Chapter 3 Using Machine Learning Insights to Target Research on the Seasonal Effects of Plants in Batch Treatment Wetlands 3.1 Introduction This chapter reviews the insights on evapotranspiration, plant effects, and saturated treatment wetlands from the in silico experiments and discusses subsequent research questions, methodology to answer them, and expected outcomes of that methodology. 3.2 Insights from machine learning on the effects of plants The machine learning work in Chapter 2 provided information about several gaps or uncertainties in the performance dynamics of treatment wetlands, particularly regarding the combined effect of plants and temperature. Plants and temperature clearly had joint impacts (Figure 2-15 and Figure 2-16) on the predicted effluent quality in the models. However, gaps in understanding the reasons and extent of these impacts were evident due to both limitations of feature inclusion and the design of treatment wetlands included in the training data. Evapotranspiration is one feature which could not be included due to insufficient data, and model interpretation of warm planted systems highlighted uncertainty in its effect. Figure 2-15 shows that plants increased predicted effluent organic matter for unsaturated systems in warmer conditions, which may have reflected evapotranspiration-caused concentration of effluent. Horizontal flow systems have previously observed this, noting effluent concentrations increasing and water levels decreasing due to summer evapotranspiration (Sultana et al., 2016), but no studies to our knowledge have quantified this effect. Despite lack of quantification, compensating for evapotranspiration has been done by multiple authors through artificial water level maintenance (Taylor et al., 2011) or mathematical 59 compensation (Sultana et al., 2016). This forced intentional omission of warm temperature results from Taylor’s work and all results from Sultana to allow any potential evapotranspiration effects to be visible in the model. Notably this also meant there were no planted warm temperature batch treatment wetlands represented in the data, as no batch studies to our knowledge allowed evapotranspiration to have an effect. This leaves considerable uncertainty around the size and type of effect that evapotranspiration has, especially for batch systems. Uncertainty around plant effects at lower temperatures was present as well. While plant status showed a large effect on the removal of organic matter at low temperatures for the overall saturated model (Figure 2-13), comparing results from a continuous flow specific model revealed a smaller effect (Figure 2-14). The subdued effect of plants in Figure 2-14 could be contributed any of operation (batch vs continuous flow), HLR, or plant species. While we know that plant species can make a difference, the impact of larger HLRs and operation on the impact of plants has not been made clear to our knowledge, with comparisons only done with significantly different HLRs (Stein et al., 2003), leaving a gap in understanding. Less represented but included in the training data were harvested treatment wetlands, which showed a slight but uncertain effect. Harvested treatment wetland data points were present in four of the publications used in the saturated training data, representing only 5% of the data points. Interpretation of harvested treatment wetlands showed either slightly worse or similar cold temperature organic removal than planted. These results matched the findings of the only publication of harvested treatment wetlands to quantify the performance impact (included in the training data) (Q. Wang et al., 2015b), but the limited data and the variability between algorithms reduces the confidence of this relationship. Other publications have 60 reported non-quantified reduced treatment after harvesting (Sultana et al., 2016) or did not mention an effect (Akratos & Tsihrintzis, 2007a; Sani et al., 2013). This reveals insufficient knowledge about the effect of harvesting for different ranges of design variables or plant species. These knowledge gaps identified through data collection and model interpretation for plants can be summarised as follows: • No information exists to our knowledge of how batch operated treatment wetlands are affected by evapotranspiration. • It is uncertain whether plants improve cold temperature treatment more in batch operated microcosms due to the batch operation, long retention times, or species selection. • Data on effluent quality for harvested treatment wetlands is limited, reducing confidence in how harvesting may impact effluent treatment at design or operational conditions other than those from the singular study. 3.3 Research questions driven by machine learning insights Allowing evapotranspiration to occur in batch-operated systems would not only concentrate effluent but also lower the water level throughout each batch. This would both extend the unsaturated period of the upper media and gradually reduce the media in contact with the wastewater as the water level is reduced. Microcosm batch experiments showed increased performance by exchanging saturation time for unsaturation time, with the best treatment at an unsaturated time of 1 day for every 2 days of saturation (Jia et al., 2010). Increased unsaturated time has also shown to increase the effective porosity in as few as three days, resulting in implications for bed hydraulics and volumetric capacity for each batch (Hua et 61 al., 2014a). Given these results, how does evapotranspiration-driven water level reduction in batch treatment wetlands affect both treatment performance and bed hydraulics? Evapotranspiration must be less than HLR to maintain wet conditions in batch operated treatment wetlands in summer. When winter comes, can plants still improve batch treatment at these higher HLRs? Evapotranspiration rates of non-aerated horizontal flow treatment wetlands have been reported to peak at 19 – 28 mm/d (Nivala et al., 2022) for Phragmites australis, depending on the wetland depth, and additionally vary (5 – 37 mm/d) depending on the wetland plant species (Milani et al., 2019). This provides a practical HLR range to determine cold temperature batch treatment wetland performance. As C. utriculata has previously been shown to have the largest influence on treatment in cold-temperature microcosm batch studies at low HLRs (3 - 7 mm/d) (W. C. Allen et al., 2002b; Taylor et al., 2011), it is a suitable species to determine the effect of these larger HLRs. By using a species with known cold temperature benefit in batch operation, we eliminate two of three unknowns (species, HLR, and operation) in the knowledge gap and ask if C. utriculata specifically can maintain its cold-temperature benefits under HLRs which prevent drying out treatment wetlands in summer. To our knowledge, C. utriculata has not been harvested in any published research and therefore represents an ideal species to widen the literature on how harvesting affects treatment. Harvesting of P. australis has been reported to worsen cold temperature performance of planted horizontal saturated wetlands at influent concentrations as high as 4000 mg/L COD at 12 mm/d (Sultana et al., 2016), as well as influent as low as 60 mg/L COD for a 7-day batch (Q. Wang et al., 2015b). As P. australis has been shown to have a less pronounced cold temperature benefit than C. utriculata (Taylor et al., 2011), will C. 62 utriculata have a larger drop in treatment when harvested, with more to lose, or have a smaller change that gives insights on the potential mechanisms behind treatment improvements? In summary, the research questions prompted by the machine learning work are: • How does evapotranspiration-driven water level reduction in batch treatment wetlands affect both treatment performance and bed hydraulics? • Can C. utriculata maintain its cold-temperature benefits under HLRs which prevent drying out treatment wetlands in summer? • 3.4 How will cold-temperature treatment performance change if C. utriculata is harvested? Methodological decisions To study these research questions, validate observations of plant effects from model interpretation, and enable future inclusion of planted batch wetland training data in warm conditions, an experiment using batch treatment wetland microcosms was designed. A HLR of 22 mm/d was chosen to match the maximum evapotranspiration observed in batch treatment wetland microcosm research (W. C. Allen et al., 2002a) while staying within the reported ranges of evapotranspiration for larger treatment wetlands (Nivala et al., 2022). C. utriculata was planted in half the microcosms to make its effect clear for both us and for training future models. Both total and soluble biochemical oxygen demand (COD) were included as measured parameters to clarify differences in temperature impact on soluble or total measurements of organic matter. While not part of the main discussion, soluble COD had a smaller but uncertain response to colder temperatures (Figure B-3) in model interpretation of saturated 63 systems. The SVM and random forest algorithms showed soluble COD to respond less to temperature, while the other algorithms were closer to parallel with total measurements of organic matter. This may be partially due to limited publications in the training data which reported both soluble and total organic matter, limiting the model’s capacity to concretely determine the relationship between them. Investigating this relationship further would also improve understanding of the potential differences between studies which treat wastewater with particulate and those with only solubilized COD sources. To fully understand the potential impact of evapotranspiration on treatment performance, mass and concentration of COD and NH3-N effluent parameters were measured. Mass provides an accurate measure of how much of each parameter is actually being removed during high evapotranspiration (Sultana et al., 2016), allowing potential effects of increased unsaturated time to be more visible. Concentration is important as a regulated parameter (BC Reg 87/2012, 2012; SOR/2012-139, 2012), with NH3-N and COD both able to be integrated into training future machine learning models. As evapotranspiration may also affect the effective porosity, water level and mass measurements were required. By monitoring the volume of water required to fill each batch microcosm to a set depth, as well as the water level at the end of the batch and the volume drained, the effective porosity at the start and end of each batch can be measured. Effective porosity has been observed to change substantially over time (Hijosa-Valsero et al., 2012; H. Liu et al., 2018), and be a valuable machine learning predictor for treatment performance (Akratos et al., 2008). Effective porosity correlates directly with clogging potential (Zhao et al., 2009), which continues to be a challenging phenomenon to predict. Clogging has similar modelling problems to treatment performance, in that prediction methods are either simple 64 rules of thumb such as HLR (Knowles et al., 2011) or require many detailed parameters for complex models (Fang et al., 2022). Effective porosity was not reported sufficiently to allow inclusion in this research’s machine learning, however clogging research may be assisted by machine learning in the future if enough publications are able to measure it. 3.5 Expected outcomes, connections, and significance This targeted research post-machine learning can reduce uncertainties around the machine learning predictions of how plants and batch operation affect treatment. This includes the effect of effluent temperature on treatment performance based on soluble or total COD (Figure B-3), operation (batch versus continuous flow), and planted, harvested, or unplanted status (Figure 2-13, Figure 2-14, Figure 2-16). By adding measurements of parameters such as water level and effective porosity, the possible mechanisms causing the effluent concentrations can be further discussed, furthering intelligent feature inclusion and analysis of interpretation efforts in modelling. This experimentation of a specific subset of the machine learning data allows furthering of design and operational recommendations using their combined insights. Machine learning work is important for identifying and providing initial quantification of design and operational effects. Unfortunately, less common design combinations, specific items like plant species, or less reported parameters such as effective porosity are currently not able to be effectively included, and require further experimentation to determine their effect. Targeted research can reveal the dynamics of these parameters to enable their inclusion in future machine learning, develop an improved understanding of specific treatment wetland configurations, and further design and operational recommendations for colder climates. 65 Chapter 4 C. utriculata Improves Treatment in Three-Day Batch Operated Treatment Wetlands in Warm, Cold, And Harvested Conditions 4.1 Introduction Horizontal subsurface constructed wetlands suffer declining treatment capabilities as water temperatures lower (Chouinard et al., 2015; Varma et al., 2021). This can be problematic for seasonally colder climates where regulations on effluent must be met year-round. Fortunately batch constructed wetland microcosms have shown that particular plant species have high levels of root zone oxidation when both air and water temperatures are at 4°C, greatly improving chemical organic matter (COD) and ammonia (NH4) removal (C. R. Allen et al., 2023; Taylor et al., 2011). In these studies, Carex utriculata has shown to have the best coldtemperature performance. A potential limitation of existing batch studies of cold temperature planted systems are their low hydraulic loading rates (4 - 7 mm/d) compared to evapotranspiration (22 mm/d) (W. C. Allen et al., 2002c; Taylor et al., 2011). In real systems, shorter batch times would be required to avoid drying out plants, and water levels would still decrease during the batch. These shorter batches may lower plant benefits during cold conditions (Taylor et al., 2011), while warm temperature water level decrease would extend the unsaturated period, potentially improving treatment (Jia et al., 2010) and preventing or delaying long-term clogging, as measured by effective porosity (Hua et al., 2014a). Increased unsaturated periods may also increase the root to shoot ratio (J. Zhang et al., 2024), altering root dynamics that impact performance even once evapotranspiration slows in colder weather. 66 Particulate sources of COD may also affect how cold temperature, and plants jointly affect treatment performance. Existing batch research has used hydrolyzed sources of COD, which are more readily biodegradable than the substantial percentages of particulate COD seen in real municipal and domestic wastewater (Gujer et al., 1999). In cold temperatures hydrolysis is reduced (Y. Huang et al., 2005), resulting in an accumulation of slowly biodegradable particulate COD (Samsó & Garcia, 2013). Plant roots can increase total suspended solids removal (Karathanasis et al., 2003), but the overall impact of plant species is not yet fully understood in these conditions. Finally, reported cold temperature benefits may conflict with the practice of harvesting plants to remove sequestered nutrients. Nutrient uptake is important for less aerobic systems such as horizontal constructed wetlands since oxygen sources for nitrification are minimal, and phosphorus uptake is important as adsorption is finite (Vymazal, 2005). To prevent sequestered nutrients from returning to the wetland through the winter and following spring by leaching out of decomposing vegetation, fall cutting and removal of the above ground plant stalk is practiced (Vymazal, 2020). Unfortunately, harvesting Phragmites australis has been shown to reduce root zone oxidation in low-strength batch microcosms (Q. Wang et al., 2015a), and increase effluent concentrations in horizontal wetlands treating high-strength dairy wastewater (Sultana et al., 2016). No studies of harvesting C. utriculata have been conducted, nor studies comparing unplanted and planted systems, which are both of interest due to previously reported high cold temperature performance. To address the current research gaps, planted (C. utriculata) and unplanted microcosms were operated in three-day batches through warm and cold conditions, treating a synthetic wastewater that included particulate COD. Treatment metrics measured included total COD 67 (tCOD), soluble COD (sCOD), NH4, PO4, pH, and redox. Warm condition operation determined how evaporative drawdown over shorter batch times alters how plants affect both removal efficiencies and the effective porosity. Cold temperature operation determined how previously reported treatment benefits of C. utriculata changed with more realistic batch times. Finally, harvesting C. utriculata quantified how harvesting impacts treatment performance for this species, and if the harvested systems still outperformed unplanted. 4.2 Materials and Methods 4.2.1 Experimental setup Three replicates each of planted (C. utriculata) and unplanted batch-operated fill-and-drain microcosms were constructed from 0.22 m2, 0.30 m deep HDPE bins. Each bin was constructed with a lower 40 mm layer of 38 mm diameter gravel and an upper 260 mm layer of 9.5 mm diameter round gravel, separated by a fiberglass mesh screen. Ratchet straps and wooden boards were used around the exterior of each bin to prevent long-term bulging of bin sides that could alter the area and depth of the bin. A 100 mm diameter sampling pipe with 100 mm spaced 6 mm holes was installed vertically in the middle of each system. For the three planted systems, C. utriculata was planted in plastic mesh pots (130 mm diameter x 100 mm depth) with six to eight stalks of mixed maturity per pot and buried in the top 100 mm of the bed at six pots per planted system. 68 Figure 4-1: Planted microcosm detail drawing. Note that unplanted microcosms are the same, minus the planted pots. All gravel and plants were sourced from previous microcosms which had been operating since June 2023 (10 months prior) with varied temperature (6°C, 20°C), influent strength (2000 mg/L tCOD, 300 mg/L TN; 500 mg/L tCOD, 75 mg/L TN), and plant (unplanted, C. utriculata) factors. All gravel and plants were mixed prior to use within the new microcosms to ensure homogenous starting conditions but were not washed to encourage a quicker startup time. 4.2.2 Synthetic wastewater composition Synthetic wastewater was prepared using dog food concentrate, urea fertilizer concentrate, and reverse osmosis water. The dog food concentrate was prepared as a modified version of the procedure developed by Kargol et al. (2023) by blending dog food (Purina Dog Chow Complete Adult Dog Food with Real Chicken) with warm tap water at 90 g/L for two 69 minutes, pouring it through a coarse sieve, swirling to remove foam, allowing the mixture to soak and settle for 12 – 24 hours, and decanting the top 50%. By adding the steps of blending and decanting, COD concentration increased to 50,000 ± 11,000 mg/L with a higher ratio of particulate COD. This also resulted in 64% less dog food used compared to Kargol’s procedure, and their noted 67% decrease in COD from dog food stored for longer than three months was not experienced in this study. Urea fertilizer concentrate was prepared by dissolving granular urea fertilizer (Nutrien Ag Solutions Evergro Special Fertilizer Urea Nitrogen 46-0-0) with DI water for a urea-N concentration of 13,800 mg/L. The final wastewater was prepared immediately before feeding by mixing the dog food and urea concentrates with reverse osmosis water to meet target concentrations (see Table 4-1). Reverse osmosis water was sourced from tanks in the enhanced forestry lab (EFL) or cold room depending on the experiment phase. Table 4-1: Synthetic wastewater composition Parameter Total COD Soluble COD Phosphate-P Total nitrogen pH 4.2.3 Value 509 ± 113 155 ± 19 3.9 ± 1.3 52.6 ± 1.7 7.0 ± 0.1 Unit mg/L mg/L mg/L mg/L Operation, sampling, and measurements Microcosms were operated in batch mode. For each batch, microcosms were filled using buckets of synthetic wastewater to a depth of 0.25 m, left for three days, drained for one hour into buckets, and refilled with new wastewater. Water level was measured at the start and end of each batch. The volume dosed and drained from each microcosm was measured by weighing the buckets before and after dosing and draining. For each microcosm, the drained 70 wastewater was stirred to ensure homogeneity, tested for soluble COD using an APURE COD-208 COD sensor (UV-Vis absorption: Xu et al., 2023), and sampled (2 x 45 ml per system). Samples were measured for total COD (closed reflux 17 colorimetric: APHA 2017: 5220-D), NH4-N (ammonia selective electrode: APHA 2017: 4500 D), NO3-N (cadmium reduction method: APHA 2017: 4500 NO3-B), PO4-P (ascorbic acid: APHA 2017: 4500 P-E), and pH (electrode: APHA 2017: 2580-B). Continuously monitoring probes within one planted and one unplanted microcosm measured bed temperature (probe method: APHA 2017: 2550-B, pH (electrode: APHA 2017: 2580-B), ORP (Pt electrode method: APHA 2017: 2580-B), and DO (membrane electrode: APHA 2017: 4500 O-G). Air temperature was measured by room sensors. The effective porosity was calculated in two ways: a fill porosity calculated from the volume of wastewater dosed to achieve a 0.25 m depth, and a drain porosity calculated from the volume drained at batch end and the depth of effluent immediately prior to draining. Figure 4-2: From left to right: microcosms in the EFL (batches 0 – 29), the cold room (batches 29 – 50), and after being harvested in the cold room (batches 50 – 70) Microcosms were first matured in a climate-controlled EFL with natural summer lighting conditions from April 30, 2024 to July 26, 2024 at the University of Northern British Columbia, Prince George, Canada. Once both plant size and effluent quality was stabilized (judged based on similar performance for three consecutive batches), microcosms were 71 moved to a cold room with artificial lighting on a 12 hour on/12 hour off cycle, which provided photosynthetically active radiation (PAR) measured at 170 μmol/m²/s using a handheld PAR sensor. Once water quality had again stabilized, plant stalks were cut 1 – 2 cm above the gravel, and monitoring continued until a new stabilization point had been reached, at which point the experiment was concluded. Table 4-2: Air and water temperature throughout the experiment, depending on the location of the microcosms. Values are listed as mean ± standard deviation (min – max) to show variability. Parameter Air temperature Water temperature EFL [°C] 21.6 ± 3.0 (16.5 – 27.1) 23.0 ± 2.5 (18.0 – 30.7)* Cold room [°C] 4.6 ± 1.1 (2.3 – 8.0) 5.5 ± 0.5 (4.3 – 8.0) *Note that water temperature in the EFL is only from a single planted system. 4.3 Results and discussion 4.3.1 Evapotranspiration To our knowledge, evapotranspiration has not been reported for C. utriculata prior to this study. By the end of the warm conditions, evapotranspiration (Figure 4-3) was within the same range as outdoor microcosms (15 – 22 mm/d) (Hijosa-Valsero et al., 2012) and pilotscale continuous flow experiments (18 – 28 mm/d maximums) (Nivala et al., 2022) which used the more commonly reported P. australis. Warm versus cold temperature planted evapotranspiration rates matched the ratio of the Jackson Meadow horizontal aerated wetland system evapotranspiration rates at similar air temperatures (Wallace & Nivala, 2005). Postharvest, differences in water loss between planted and unplanted microcosms were not statistically significant. 72 Figure 4-3: Evapotranspiration (for planted units) and evaporation (for unplanted units) throughout each batch of the experimental period. The first dotted line marks the start of cold conditions, while the second marks when plants were harvested. 4.3.2 Chemical Oxygen Demand In warm conditions mass removal of COD was identical for planted and unplanted microcosms COD (96±2%) (Figure 4-4). Longer (10 – 20 day batches) constant water level batch studies report poorer mass removal (83%) in C. utriculata microcosms compared to unplanted (88%) in warm conditions (Taylor et al., 2011). The cause of this is uncertain. Plants can add carbon to treatment wetlands through exudates and detritus, potentially worsening COD removal. Monitoring a free water surface treatment wetland for 56 days showed that 5–8% of the total Typha biomass was leached as DOC, and 30–50% exited as particulate organic carbon or was microbially respired (Pinney et al., 2000). In planted microcosms organic carbon exudates are estimated to be at least 90 mg C/m2/d for zerocarbon influent, and especially substantial for C. utriculata in the summer compared to other species (C. R. Allen et al., 2023). If this rate is consistent, the mass addition of carbon 73 exudates would have been less, just as radial oxygen loss has been expected to also have decreased effects over shorter batch times (Taylor et al., 2011). However, Taylor argued the organic carbon in the wastewater is overwhelming larger than plant additions even at higher batch durations and did not think plant carbon addition could be the cause. Regardless of the reason, the lengthened unsaturated time in our planted microcosms may have increased COD removal. In three-day batch microcosm studies, two days of contact time followed by a whole day of draining had 97% COD removal, while three days of contact time with minimal drain time had 92% removal (Jia et al., 2010). The extended unsaturated oxidation period from the 20 – 83% water level fluctuations in our study may have similarly improved planted COD removal. Figure 4-4: NH3-N, COD, and PO4-P effluent concentrations (left) and mass (right) for planted (red) and unplanted (blue) microcosms. The first line marks the start of cold room conditions, and the second line marks the time at which C. utriculata was harvested. 74 In cold conditions, C. utriculata increased the removal of particulate COD (Figure 4-4). Cold conditions caused an immediate increase in soluble COD for all microcosms, followed by a slow but larger increase in particulate COD. This particulate COD increase was less in planted microcosms, indicating additional capture or hydrolysis. Other works show that while COD removal is not influenced by particulate in warm conditions (Caselles-Osorio & García, 2006), it is limited by the rate of hydrolysis in cold conditions (Samsó & Garcia, 2013). The subsequent accumulation of particulate produces acetic acid (Y. Huang et al., 2005). Here we saw particulate increase in unplanted, and pH of unplanted systems also declining from 7.25 to 7 (Figure 4-5). Planted microcosm pH stayed constant, but also started at a lower pH of 6.75. Figure 4-5: pH in the effluent of each batch. The first line marks the start of cold room conditions, and the second line marks the time at which C. utriculata was harvested. The COD removal shown by C. utriculata in cold conditions was not as high as comparable studies with longer batches. In the relatively stable four batches before harvesting, COD mass removal was 86±5% for planted, compared to 73±9% for unplanted. Similar studies with 20day batches report 96 – 99% for C. utriculata and 58 – 67% for unplanted (Stein & Hook, 2005; Taylor et al., 2011). In planted treatment wetlands, oxygen transported through the root 75 aerenchyma diffuses through roots as radial oxygen loss (Armstrong et al., 2000), which varies substantially by species (Lai et al., 2012). Assuming a relatively consistent flow of oxygen, shorter batch times, as with our microcosms, have been thought to result in a lower mass of oxygen per batch for heterotrophic usage (Taylor et al., 2011). Maturation stage is also a possible contributor, as microcosms in the present study matured for three months prior to the cold conditions, compared to six months or more in other studies. In our machine learning models, age showed improved organic matter removal until 4 months to a year, depending on the algorithm (Figure 2-18). Microcosm research of Juncus effusus classified maturity as post-80 days, with significant differences in COD, NH3-N, and DO concentrations attributed to maturity (Hu et al., 2023). Mature specimens were transplanted but may still have impacted results due to less time for extensive root development. Harvested microcosms retained their mass and concentration advantage over unplanted microcosms. Post-harvesting, mass COD removal in harvested microcosms immediately decreased from 83% to 79%. These are proportionally similar results to a 7-day batch experiment with P. australis using 59 mg/L COD influent which reported an immediate concentration removal change from 56% to 50% upon harvest (Q. Wang et al., 2015a). Despite this change, COD mass removal remained worse in unplanted microcosms, at 72% during the same batch. Radial oxygen loss has been shown to continue from airflow through cut dead plant culms (Tanaka et al., 2007). Post-harvesting, the planted microcosms continued to improve while unplanted performance worsened. As a result, the experiment concluded with a larger difference between harvested and unplanted microcosms than existed before harvesting. 76 Evapotranspiration reduction in cold and then harvested conditions tempered the effluent concentration change for planted microcosms. When microcosms were moved to cold conditions, evapotranspiration immediately lowered from 15 – 20 mm/d to less than 5 mm/d, doubling effluent volume. As a result, the COD effluent concentration increased by only 25%, as opposed to the 400% increase in COD mass. Upon harvesting, evapotranspiration reduced again to 0.7 mm/d, causing the COD effluent concentration to increase by 30% instead of 50%. 4.3.3 Ammonia Planted microcosms had a consistently higher mass removal rate of ammonia through warm conditions (Figure 4-4), with parallel improvement in removal between planted and unplanted. Throughout this whole warm period the ammonia mass removal of the planted microcosms was consistently 50% higher. This would align with constant nutrient uptake observed during plant growth periods (Bender et al., 2015), but the literature varies. Research with Juncus effusus indicates that nitrogen removal increases once plants are mature, and uptake itself is only 2 – 4% of nitrogen removal in treatment wetlands (Hu et al., 2023). Lengthened rest periods, as observed in our microcosms, can also greatly enhance nitrification (J. Zhang et al., 2024). Ammonia removal was not correlated with the rise in evapotranspiration from 8 to 20 mm/d, however effluent mass was more concentrated, resulting in 25% lower effluent ammonia concentrations in the planted microcosms. Planted microcosms retained their advantage throughout the cold period, but the advantage was less than studies with longer batch times. Cold temperatures decrease plant metabolism, and while this results in higher root zone oxidation for nitrification (Mburu et al., 2012), it has been suggested to reduce plant uptake to a greater degree (Yates et al., 2016). Despite 77 this, our cold planted ammonia removal was better than cold unplanted microcosms and similar to warm unplanted microcosms. Overall, the ammonia removal rates were lower than similar studies, which we attribute to our lower batch time and higher organic loading. Nitrate was zero throughout the experiment for planted and unplanted, as the conditions for denitrification were better than those for nitrification. We suspect studies with high removal rates within three days of a longer batch are a result of rapid sorption of ammonia onto the biofilm (Tanner et al., 1999). Plant uptake and nitrification of this sorbed ammonia would then continue throughout the rest of the batch, taking advantage of the whole batch time. Other studies with higher removal rates had lower organic loadings. Ammonia removal decreases as COD influent increases, since this increases competition for oxidation sites (Riley et al., 2005). A less explained comparison is with horizontal flow wetlands using a 3.4 day HRT planted with Carex aquatilis, achieving 98% NH3 removal at 5 - 10°C (Yates et al., 2016). As their unplanted removal (2%) is similar to the present study, the improved planted performance suggests an interaction between their shallower depth and plant roots to potentially be of benefit, which would align with the increased evapotranspiration observed in shallower horizontal flow treatment wetlands (Nivala et al., 2022). Harvesting lowered the advantage of planted systems from 40% less ammonia mass than unplanted to only 12% less. We attribute this gradual decline to the saturation of adsorption sites (Vymazal, 2007) after plant uptake and/or oxidation decreased. The remaining benefit may be due to oxidation provided through cut plant culms (Tanaka et al., 2007) or continued accumulation of nitrogen into the plant biomass (X. Huang et al., 2020). Continued 78 accumulation has been observed for senesced Phalaris arundinacea but has not yet been proven for cut plants nor C. utriculata. 4.3.4 Phosphorus Planted microcosms had consistently low effluent phosphate concentrations, reaching 0 mg/L by batch 22 – the same time that unplanted phosphate removal started to stabilize. Once zero phosphate in the effluent had been achieved in the microcosms it was maintained regardless of temperature or plant harvesting. At stability, planted microcosms reduced phosphate by an additional 30 mg/m2/d compared to unplanted (Figure 4-4). Phosphate removal is difficult to compare to other studies as removal rates vary significantly and mesocosms significantly outperform other treatment wetlands (Ury et al., 2023). That said, the rate of 2.2 g PO4-P/m2/year dosed to our microcosms is comparable to the 2.0 g TP/m2/y Ury found across all wetlands which were phosphorus sinks, and less than the 4 g P/m2/y that can be removed through harvesting biomass in full scale horizontal subsurface flow treatment wetlands (Okada & Vymazal, 2023). As unsaturated periods have been shown to enhance inorganic phosphorus solubilization leading to easier absorption by plants (J. Zhang et al., 2024), our planted microcosms may have been able to treat a higher loading of P while maintaining low effluent concentrations. The lack of clear temperature effect on phosphate removal aligns with some literature (Sani et al., 2013), while contrasting with findings which found reductions in removal at lower temperatures (Akratos & Tsihrintzis, 2007b). The reason for our continued 0 mg/L phosphate in planted microcosms in cold temperatures and post harvest is uncertain. Adsorptive capacity may have remained due to previous plant uptake, leaving sites for phosphate adsorption post-harvest that would eventually fill and reduce removal unless plant growth 79 restarted. Continued uptake of phosphates into the root biomass has been observed in senesced (but not harvested) wetland plants throughout winter (X. Huang et al., 2020). Plant roots can also increase adsorption sites; iron plaque formed on roots has shown to be a significant contributor to phosphorus removal (Y. Liu et al., 2024) in floating treatment wetlands. 4.3.5 Porosity In warm conditions, effective porosities based on the volume to fill each batch (‘fill’ porosity) were higher in planted microcosms (Figure 4-6). Interestingly, the porosity of the planted microcosm’s lower drained portion (‘drain’ porosity) was 7 – 10% lower than the overall porosity of the entire filled microcosm. In contrast, unplanted microcosms had only a 1% difference between the ‘fill’ and ‘drain’ porosities. While not observed in other fill-anddrain treatment wetlands to our knowledge, periodic unsaturated rest periods of continuous flow microcosms have increased effective porosity by up to 6% over 10 days due to the inflow of oxygen that accelerates organic matter decomposition (Hua et al., 2014b; H. Wang et al., 2021). The reduced porosity observed in lower portions of the microcosm is therefore attributed to lower oxygen availability at depth, with the higher porosity in upper portions of the planted microcosms attributed to lengthened unsaturated periods caused by evapotranspiration. The effective ‘fill’ porosity of the planted and unplanted microcosms diverged slightly throughout the maturation period, with planted microcosms increasing and unplanted decreasing for a 2 – 3% difference by the end of the warm period. In the fill-and-drain batch configuration demonstrated in this study, resting periods would occur naturally from evapotranspiration each warm season, potentially preventing long term biological-based 80 clogging effects. Measurements of effective porosity in a long-term study would be of benefit to confirm any trends shown in our short-term study. Figure 4-6: Effective porosity of each batch, based on the volume to fill the drained microcosms to a set water level. The first line marks the start of cold room conditions, and the second line marks the time at which C. utriculata was harvested. In cold conditions effective porosity decreased substantially, and harvesting did not have a discernable impact. Thicker biofilm growth occurs in attached-growth and planted systems in lower temperatures (Morcillo & Manzanera, 2021). Lower temperatures decrease the biodegradation rate, resulting in reactants diffusing further through biofilm before being consumed and allowing a thicker biofilm to accumulate (Ahmed & Delatolla, 2021). Build up of non-degraded particulate matter can also contribute to clogging (Zhao et al., 2009), and we observed a larger fraction of particulate COD in the effluent in cold conditions. The cause of continued decline of planted microcosm porosity in the final four batches of planted microcosms (as opposed to the unplanted) is uncertain and deserves future study. This reduction in porosity between warm and cold conditions represents a substantial 17% (unplanted) to 29% (planted) reduction in working volume and mass loading over a period of 60 days. For batch operated microcosms used for experiments, this shows the importance of volumetric inflow and outflow measurements for mass removal calculations. In real world batch operated treatment wetlands, this could result in a critical reduction in flow capacity 81 unless batch length is dynamically adjusted. In the more commonly scaled-up horizontal flow treatment wetlands, a similar reduction in porosity would reduce the hydraulic retention time at a time where treatment is already lower due to temperature, potentially further worsening treatment. This may contribute to the clogging issues documented in some vertical flow wetlands during cold spells (Langergraber et al., 2009), despite their unsaturated periods. 4.4 Conclusion The findings of this study show that batch treatment wetland microcosms planted with C. utriculata have consistently better performance than unplanted equivalents in warm, cold, and cold harvested conditions, apart from COD effluent concentrations at warm temperatures. C. utriculata microcosms had higher effective porosities than unplanted, which appear to be driven by evapotranspiration caused unsaturated periods. This potential for longterm clogging resistance should be further investigated with long term studies. Warm temperature mass COD removal of planted microcosms was identical to unplanted, unlike similar studies which reported poorer performance of planted microcosms in warm conditions. In cold conditions C. utriculata continued to improve treatment compared to unplanted, although the shorter batch time was associated with a lower benefit than previously reported. Harvesting C. utriculata resulted in a 15% increase in effluent total COD mass, and a 20% increase in NH3-N mass, but continued complete removal of phosphates. Despite these reductions, harvested systems still provided better treatment for all parameters when compared to unplanted. Negative effects on effluent concentrations were minimized due to reductions in evapotranspiration during cold temperatures and after harvesting causing a higher effluent volume. 82 Chapter 5 Conclusion 5.1 Background Machine learning and microcosm experiments in this thesis have shown that the dynamics of how both warm and cold temperatures impact treatment wetlands are highly dependant on their design and operation variables (Figure 2-12). Machine learning models trained with a wide range of published data have proven to be relatively effective at accounting for these many permutations and providing a global understanding of feature effect through interpretation tools. Models are still limited by the extent of data available, and variables insufficiently reported to include and less common design permutations led to uncertainty around more specific effects and feature interactions. These areas of uncertainty identified through the machine learning process can and should be followed up with physical experiments to effectively fill in gaps in understanding. While multiple such uncertainties were identified, evapotranspiration, operation, and plant harvesting were chosen for follow up. 5.2 Joint results of machine learning and microcosm experiments In Chapter 2, machine learning interpretation of temperature predicted that planted continuous flow saturated treatment wetlands would have a 7 – 12 mg/L lower organic matter concentration than unplanted, while the general saturated model influenced by low HLR batch microcosms showed a larger 18 – 30 mg/L benefit. Chapter 4’s physical experiments with C. utriculata revealed a 30 mg/L benefit of planted microcosms even at 22 mm/d, matching the higher end of the general saturated model effect. While not as impactful as 3 mm/d microcosms using C. utriculata (Taylor et al., 2011), it shows that a significant 83 percentage of the plants benefit is retained at higher HLRs, and inclusion of plant species will be important in future machine learning models. Adding Chapter 4’s microcosm data to the machine learning training data (Figure 5-1), results in between the original saturated model and the continuous flow-only model were predicted. Figure 5-2 shows the model assumes these plant benefits to be applicable to both fill-and-drain and permanently saturated treatment wetlands, since PDPs include all possible combinations and interactions, and saturated treatment wetlands would have poorer performance if the model understood plants to have less benefit for permanently saturated treatment wetlands. Figure 5-1: PDP curves of the effect of plants dependant on temperature AFTER Chapter 4 microcosm data was added. Results are in between the original saturated model and the continuous flow only model. Harvesting shows a larger effect than previously predicted. Warm temperature planted batch treatment wetlands were missing from Chapter 2’s data collection due to intentional exclusion treatment wetlands with artificial controls on evapotranspiration. This gap was filled by the microcosm experiments in Chapter 4. While the evapotranspiration water level reduction did not result in a definite improvement in NH3N or PO4-P removal compared to other studies, there was increased COD mass removal compared to longer 20-day batch studies using C. utriculata. In addition, planted microcosms observed a lower effective porosity attributed to the extended unsaturated time from evapotranspiration. A revised model including the new microcosm data did not substantially change the effect of operation on the average predicted organic matter effluent (Figure 5-2, as 84 compared to Figure B-4). No algorithms predicted effluent concentrations to be higher for planted fill-and-drain treatment wetlands in summer despite the new microcosm data. Figure 5-2: PDP curve for the effect of temperature depending on the whether the treatment wetland is fill-and-drain (batch) operated or permanently saturated. Model shows results AFTER microcosm data added. Particulate COD was affected by temperature to a much larger degree than soluble COD in the Chapter 4 microcosms, supporting the machine learning findings. Random forest and SVM PDPs (Figure B-3) showed soluble COD to be less effected, despite the obvious errors of soluble COD being larger than total in the PDP. With added data from the new microcosms, the modelled relationships changed to those of Figure 5-3. While sCOD is still higher than COD for warm temperatures in the random forest model, the total COD is now shown to more affected by cold temperatures than sCOD for three of the four models, better matching the experimental results. 85 Figure 5-3: PDP curves AFTER microcosm data added for the effect of temperature, depending on the organic matter parameter for saturated systems. COD is more consistently greater than sCOD in the revised plots, especially at colder temperatures. Harvesting the C. utriculata in Chapter 4’s microcosms caused an immediate increase in effluent particulate COD while soluble COD effluent remained constant. These results indicate that plant stalks were increasing the rate of hydrolyzation while having little impact on the degradation of soluble COD, suggesting the impact of plant species on soluble and total COD removal in cold temperatures should be different. As a three-way interaction, this cannot easily be visualised through current model interpretation techniques. 5.3 Cold climate design and operation When considering design in cold climate locations, the influent wastewater temperature is critical as it is heavily responsible for the resultant wastewater treatment process temperature (Figure 2-2). Once the type of treatment wetland is determined, it can be considered together with the HLR and air temperature data to model the effluent temperature (Figure 2-3), representing the coldest water in the treatment wetland (Figure 2-17). While insulation is known to be effective in other publications, quantification of this effect was not achieved in this work. Well designed forced aeration and alternating treatment wetlands are do not appear to be affected by temperature at all in the design and temperature ranges currently explored (Nivala 86 et al., 2019), but do have additional power and operational requirements which are outside the scope of this thesis. Neither design was studied in this work due to less data availability. For the other main treatment wetland types, unsaturated vertical flow treatment wetlands are less temperature sensitive than more saturated treatment wetlands such as continuous flow or fill-and-drain (Figure 2-12). Multiple treatment wetland design variables and operational methods were shown by machine learning to alter low and high temperature treatment rates similarly. A lack of clear modelled interaction between effluent temperature and HLR, media diameter, depth, or age indicates that the PDP relationships for these features apply across all temperatures. As such, for saturated treatment wetlands PDPs indicate that a low HLR (Figure 2-9), small media diameter (Figure 2-10), shallow depth (Figure 2-11), and old (Figure 2-18) treatment wetland would perform best. For unsaturated treatment wetlands, a smaller media diameter, thicker media layer, and shorter interval between loadings (Figure B-2) increase treatment performance. Features shown to affect warm and cold treatment wetlands differently included plants, operation, and particulate matter in the wastewater. For batch treatment wetlands, planting appropriate species (C. utriculata being best (Taylor et al., 2011)) can substantially improve cold temperature treatment up to at least 22 L/m2/d (Figure 4-4), and machine learning suggests these effects are applicable to continuous operated treatment wetlands as well (Figure 5-2). Plant benefits are slightly reduced if the plant stalks are cut and removed, but harvested treatment wetlands still hold a strong advantage to unplanted treatment wetlands (Figure 4-4). The fraction of particulate in wastewater is important as soluble COD is more 87 readily degradable than particulate in colder treatment wetlands (Figure 5-3), and increased solubilization appears to be a portion of the benefit provided by plants (Figure 4-4). 5.4 Limitations and recommendations for future work Careful interpretation of machine learning algorithms showed high potential for compiling and summarizing the relationships discovered by individual studies, but only if checked against literature. Unlike mechanistic models, machine learning has no concept of causality, and results in this thesis and future works must be used with caution and an understanding of what should theoretically occur in treatment wetlands. Interpretation of multiple algorithms together and the use of ICE curves helped to show the potential variability of relationships mapped by models and are recommended for future work. As treatment wetlands are high dimensional systems, having interpretation tools for only one or two features limited results. As data availability increases in the future, we recommend splitting training data into further subsets of treatment wetland types, as done here for unsaturated and saturated treatment wetlands, to allow for more specific interpretation plots. Data division and modelling must be done with an understanding of imbalance, correlation, or cardinality issues in the source data, or misleading results will occur. As more data-dense areas of interpretation plots tended to show increased agreement with established theory, the ability to learn from models will increase and improve as more data is available and added. We recommend that future treatment wetland studies should consider the features used here for machine learning as the minimum reporting requirements. This will enable machine learning models to add data from these studies and continually improve as research progresses. It will also provide sufficient information for normal comparison with other 88 studies, as parameters such as influent and effluent temperature have been proven here to be important but were not readily available in all publications. Theoretically important features not able to be added due to data availability issues should also be considered, such as snow presence and depth over the winter, organic nitrogen influent concentration, time between fill-and-drain batches, evapotranspiration, the influent distribution system details. Once sufficient treatment wetlands report these parameters, they could be added to future machine learning models with data imputation methods used to fill in the missing values in past publications without this data. This will allow for more accurate and comprehensive predictions and interpretation plots. It should be emphasized that while this project included only data from published treatment wetland studies, data from non-academic sources could easily be added if of high quality and containing sufficient features. Machine learning is a wide and evolving field, with many different possible approaches that could be explored to follow up to the in silico experiments in this project. A significant limitation to the approach in this project is the lack of accounting for within-wetland autocorrelation. Real wetlands are not necessarily steady state, dynamically responding to influent composition changes, long term development of microbial and plant communities, and seasonal temperature patterns. This was partially accounted for with the ‘age’ feature, but past states were not accounted for, including parameters important for temperature and plant interactions such as season, which posed difficulties in our smaller dataset due to high correlations with effluent temperature. As data availability increases over time and machine learning methods continue to develop, more options will emerge. Future modelling efforts could experiment with adding features such as season, past influent and temperature conditions, or even month. Methods could also include recurrent neural networks such as the 89 long short-term memory (LSTM) models that have been used for individual constructed wetlands (Yang et al., 2023), or methods as of yet unexplored for treatment wetlands such as temporal fusion transformers (Lim et al., 2020) which can deal with static and temporal variables while still being interpretable. Given the importance of interpretation demonstrated in this work, we recommend focusing efforts on interpretable modelling methods. Both in silico and ex situ experiments were limited by the microcosms included in both, whose smaller scales and/or controlled indoor conditions for often short periods can impact the relevance of findings. Treatment wetlands, though designed, are also natural, biological systems that can perform differently at larger scales and natural conditions over multiple years or decades (Brisson et al., 2024). These limitations are particularly relevant for the temperature and season focused work of this thesis, and data from pilot or full-scale treatment wetlands over multiple seasons or years can provide valuable improvement and follow up to all aspects of this work. Multi-year treatment wetland experiments through seasonally frozen conditions are recommended to follow up to our work with batch microcosms. Realistic weather patterns are difficult to simulate in indoor conditions, and multi-year studies are difficult for a master’s length project. While our microcosms experienced the same air and water temperatures, air and water temperature are rarely the same in outdoor conditions (Figure 2-17). As radial oxygen loss is a function of air temperature (Lai et al., 2012) while microbial communities are more affected by water temperatures, this is a limitation of our work. Cold climates also have long periods of frost and snow, unexplored in this work, which may alter passage of air through cut plant stalks, carbon release of decomposing tissue, and how C. utriculata impacts treatment. These cold periods are followed by warmer conditions in 90 spring; however, plant growth is immature at this point, again causing potential changes in how planted treatment wetlands behave. We also demonstrated that the effective porosity of batch fill-and-drain wetlands is affected by evapotranspiration water level reduction and temperature, but how these influence the effective porosity over multiple season changes and years of plant growth, harvesting, and frost damage is unknown. Experiments throughout multiple natural seasonal changes would extend our understanding of batch treatment wetland dynamics and better allow for temporal modelling that accounts for seasonal shifts. Finally, while the machine learning and microcosms of this project assume a monoculture, polycultures and changes in the dominant species can arise from ecological pressure over the long term (Zheng et al., 2020). As cold climate treatment performance has been shown to change with species (Taylor et al., 2011), this could pose issues for designs that rely on the performance of a particular species such as C. utriculata. How treatment wetland design and operation relate to the long term stability of a particular species or mixture of species is yet to be fully explored, but could potentially draw on knowledge from natural analogous wetlands (MacKenzie & Moran, 2004), and involve intentional planting of polycultures with different growth forms (Luo et al., 2023). As polycultures become more studied, inclusion strategies will also become necessary for machine learning. Plant species were already unable to be included in our model to due high cardinality, therefore strategies to reduce this and allow inclusion should be explored such as grouping similar performing species or polycultures. 91 References Agnieszka, K., & Anna, B. (2013). Cold season temperature pattern over a Horizontal subsurface flow constructed wetland—Application of thermography in environmental research. Fresenius Environmental Bulletin, 22, 3694–3699. Aguirre, P., Ojeda, E., García, J., BarragÁn, J., & Mujeriego, R. (2005). Effect of Water Depth on the Removal of Organic Matter in Horizontal Subsurface Flow Constructed Wetlands. Journal of Environmental Science and Health, Part A, 40(6–7), 1457– 1466. https://doi.org/10.1081/ESE-200055886 Ahmed, W., & Delatolla, R. (2021). Biofilm and microbiome response of attached growth nitrification systems across incremental decreases to low temperatures. Journal of Water Process Engineering, 39, 101730. https://doi.org/10.1016/j.jwpe.2020.101730 Akratos, C. S., Papaspyros, J. N. E., & Tsihrintzis, V. A. (2008). An artificial neural network model and design equations for BOD and COD removal prediction in horizontal subsurface flow constructed wetlands. Chemical Engineering Journal, 143(1), 96– 110. https://doi.org/10.1016/j.cej.2007.12.029 Akratos, C. S., & Tsihrintzis, V. A. (2007a). Effect of temperature, HRT, vegetation and porous media on removal efficiency of pilot-scale horizontal subsurface flow constructed wetlands. Ecological Engineering, 29(2), 173–191. https://doi.org/10.1016/j.ecoleng.2006.06.013 Akratos, C. S., & Tsihrintzis, V. A. (2007b). Effect of temperature, HRT, vegetation and porous media on removal efficiency of pilot-scale horizontal subsurface flow 92 constructed wetlands. Ecological Engineering, 29(2), 173–191. https://doi.org/10.1016/j.ecoleng.2006.06.013 Allen, C. R., Burr, M. D., Camper, A. K., Moss, J. J., & Stein, O. R. (2023). Seasonality, C:N ratio and plant species influence on denitrification and plant nitrogen uptake in treatment wetlands. Ecological Engineering, 191, 106946. https://doi.org/10.1016/j.ecoleng.2023.106946 Allen, W. C., Hook, P. B., Biederman, J. A., & Stein, O. R. (2002a). Temperature and Wetland Plant Species Effects on Wastewater Treatment and Root Zone Oxidation. Journal of Environmental Quality, 31(3), 1010–1016. https://doi.org/10.2134/jeq2002.1010 Allen, W. C., Hook, P. B., Biederman, J. A., & Stein, O. R. (2002b). Temperature and Wetland Plant Species Effects on Wastewater Treatment and Root Zone Oxidation. Journal of Environmental Quality, 31(3), 1010–1016. https://doi.org/10.2134/jeq2002.1010 Allen, W. C., Hook, P. B., Biederman, J. A., & Stein, O. R. (2002c). Temperature and Wetland Plant Species Effects on Wastewater Treatment and Root Zone Oxidation. Journal of Environmental Quality, 31(3), 1010–1016. https://doi.org/10.2134/jeq2002.1010 Arias, C. A., Istenic, D., Stein, O., Zhai, X., Kilian, R., Vera-Puerto, I., & Brix, H. (2022). Effects of effluent recycle on treatment performance in a vertical flow constructed 93 wetland. Ecological Engineering, 180, 106675. https://doi.org/10.1016/j.ecoleng.2022.106675 Armstrong, W., Cousins, D., Armstrong, J., Turner, D. W., & Beckett, P. M. (2000). Oxygen Distribution in Wetland Plant Roots and Permeability Barriers to Gas-exchange with the Rhizosphere: A Microelectrode and Modelling Study with Phragmites australis. Annals of Botany, 86(3), 687–703. https://doi.org/10.1006/anbo.2000.1236 Awad, M., & Khanna, R. (2015). Support Vector Regression. In M. Awad & R. Khanna (Eds.), Efficient Learning Machines: Theories, Concepts, and Applications for Engineers and System Designers (pp. 67–80). Apress. https://doi.org/10.1007/978-14302-5990-9_4 Ayotte, S. H., Wallace, S. J., Allen, C. R., Weber, K. P., Stein, O. R., & Lauchnor, E. G. (2024). Microbial community dynamics in a two-stage treatment wetland: Insights from treating seasonal ski resort wastewater. Bioresource Technology Reports, 27, 101885. https://doi.org/10.1016/j.biteb.2024.101885 Beebe, D. A., Castle, J. W., Molz, F. J., & Rodgers, J. H. (2014). Effects of evapotranspiration on treatment performance in constructed wetlands: Experimental studies and modeling. Ecological Engineering, 71, 394–400. https://doi.org/10.1016/j.ecoleng.2014.07.052 Bender, R. R., Haegele, J. W., & Below, F. E. (2015). Nutrient Uptake, Partitioning, and Remobilization in Modern Soybean Varieties. Agronomy Journal, 107(2), 563–573. https://doi.org/10.2134/agronj14.0435 94 Bonaccorso, G. (2018). Machine Learning Algorithms: Popular algorithms for data science and machine learning, 2nd Edition. Packt Publishing Ltd. Brisson, J., Carvalho, P., Stein, O., Weber, K., Brix, H., Zhao, Y., & Zurita, F. (2024). Smallscale experiments: Using mesocosms and microcosms for testing hypotheses in treatment wetland research. Ecological Engineering, 208, 107378. https://doi.org/10.1016/j.ecoleng.2024.107378 Burgos, V., Araya, F., Reyes-Contreras, C., Vera, I., & Vidal, G. (2017). Performance of ornamental plants in mesocosm subsurface constructed wetlands under different organic sewage loading. Ecological Engineering, 99, 246–255. https://doi.org/10.1016/j.ecoleng.2016.11.058 Caselles-Osorio, A., & García, J. (2006). Performance of experimental horizontal subsurface flow constructed wetlands fed with dissolved or particulate organic matter. Water Research, 40(19), 3603–3611. https://doi.org/10.1016/j.watres.2006.05.038 Chen, T., & Guestrin, C. (2016). XGBoost: A Scalable Tree Boosting System. Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 785–794. https://doi.org/10.1145/2939672.2939785 Chouinard, A., Anderson, B. C., Wootton, B. C., & Huang, J. J. (2015). Comparative study of cold-climate constructed wetland technology in Canada and northern China for water resource protection. Environmental Reviews, 23(4), 367–381. 95 Cui, D., Liang, S., & Wang, D. (2021). Observed and projected changes in global climate zones based on Köppen climate classification. WIREs Climate Change, 12(3), e701. https://doi.org/10.1002/wcc.701 Ding, J., Jia, Y., Zhao, C., Bo, W., Xu, X., Lv, R., Zhou, G., Kong, Q., Du, Y., Xu, F., & Wang, Q. (2021). Microbial abundance and community in constructed wetlands planted with Phragmites australis and Typha orientalis in winter. International Journal of Phytoremediation, 23(14), 1476–1485. https://doi.org/10.1080/15226514.2021.1907737 Dotro, G., Langergraber, G., Molle, P., Nivala, J., Puigagut, J., & Stein, O. (2017). Biological Wastewater Treatment Series, Volume 7: Treatment Wetlands. Fang, Y., Kong, L., Zhang, P., Zhang, L., Zhao, H., Xiang, X., Cheng, S., Zhang, H., Ju, F., & Li, L. (2022). Fifteen-year analysis of constructed wetland clogging: A critical review. Journal of Cleaner Production, 365, 132755. https://doi.org/10.1016/j.jclepro.2022.132755 Farrar, D. E., & Glauber, R. R. (1967). Multicollinearity in Regression Analysis: The Problem Revisited. The Review of Economics and Statistics, 49(1), 92–107. https://doi.org/10.2307/1937887 Fisher, A., Rudin, C., & Dominici, F. (2019). All Models are Wrong, but Many are Useful: Learning a Variable’s Importance by Studying an Entire Class of Prediction Models Simultaneously (No. arXiv:1801.01489). arXiv. https://doi.org/10.48550/arXiv.1801.01489 96 Foladori, P., Ruaben, J., & Ortigara, A. R. C. (2013). Recirculation or artificial aeration in vertical flow constructed wetlands: A comparative study for treating high load wastewater. Bioresource Technology, 149, 398–405. https://doi.org/10.1016/j.biortech.2013.09.099 Friedman, J. H. (2001). Greedy function approximation: A gradient boosting machine. The Annals of Statistics, 29(5), 1189–1232. https://doi.org/10.1214/aos/1013203451 Friedman, J. H., & Popescu, B. E. (2008). Predictive learning via rule ensembles. The Annals of Applied Statistics, 2(3). https://doi.org/10.1214/07-AOAS148 Garcıá , J., Ojeda, E., Sales, E., Chico, F., Pıŕ iz, T., Aguirre, P., & Mujeriego, R. (2003). Spatial variations of temperature, redox potential, and contaminants in horizontal flow reed beds. Ecological Engineering, 21(2), 129–142. https://doi.org/10.1016/j.ecoleng.2003.10.001 Grebenshchykova, Z., Brisson, J., Chazarenc, F., & Comeau, Y. (2020). Two-year performance of single-stage vertical flow treatment wetlands planted with willows under cold-climate conditions. Ecological Engineering, 153, 105912. https://doi.org/10.1016/j.ecoleng.2020.105912 Grebenshchykova, Z., Forquet, N., Brisson, J., Lachapelle-T, X., Comeau, Y., & Chazarenc, F. (2023). Thermal properties of treatment wetlands operated under freezing conditions. Water Science and Technology: A Journal of the International Association on Water Pollution Research, 88(8), 2054–2067. https://doi.org/10.2166/wst.2023.320 97 Greenwell, B. (2024). fastshap: Fast Approximate Shapley Values. https://github.com/bgreenwell/fastshap Gujer, W., Henze, M., Mino, T., & Loosdrecht, M. van. (1999). Activated sludge model No. 3. Water Science and Technology, 39(1), 183–193. https://doi.org/10.1016/S02731223(98)00785-9 Gunasekaran, A., Mistry, P., & Chen, M. (2024). Which Explanation Should be Selected: A Method Agnostic Model Class Reliance Explanation for Model and Explanation Multiplicity. SN Computer Science, 5(5), 503. https://doi.org/10.1007/s42979-02402810-8 Hassija, V., Chamola, V., Mahapatra, A., Singal, A., Goel, D., Huang, K., Scardapane, S., Spinelli, I., Mahmud, M., & Hussain, A. (2024). Interpreting Black-Box Models: A Review on Explainable Artificial Intelligence. Cognitive Computation, 16(1), 45–74. https://doi.org/10.1007/s12559-023-10179-8 He, Q., & Mankin, K. (2007). Performance Variations Of Cod And Nitrogen Removal By Vegetated Submerged Bed Wetlands. Journal of the American Water Resources Association, 01176. https://onlinelibrary.wiley.com/doi/abs/10.1111/j.17521688.2002.tb04373.x Headley, T. R., Herity, E., & Davison, L. (2005). Treatment at different depths and vertical mixing within a 1-m deep horizontal subsurface-flow wetland. Ecological Engineering, 25(5), 567–582. https://doi.org/10.1016/j.ecoleng.2005.07.012 98 Hijosa-Valsero, M., Sidrach-Cardona, R., & Bécares, E. (2012). Comparison of interannual removal variation of various constructed wetland types. Science of The Total Environment, 430, 174–183. https://doi.org/10.1016/j.scitotenv.2012.04.072 Hu, X., Yue, J., Yao, D., Zhang, X., Li, Y., Hu, Z., Liang, S., Wu, H., Xie, H., & Zhang, J. (2023). Plant development alters the nitrogen cycle in subsurface flow constructed wetlands: Implications to the strategies for intensified treatment performance. Water Research, 246, 120750. https://doi.org/10.1016/j.watres.2023.120750 Hua, G., Zeng, Y., Zhao, Z., Cheng, K., & Chen, G. (2014a). Applying a resting operation to alleviate bioclogging in vertical flow constructed wetlands: An experimental lab evaluation. Journal of Environmental Management, 136, 47–53. https://doi.org/10.1016/j.jenvman.2014.01.030 Hua, G., Zeng, Y., Zhao, Z., Cheng, K., & Chen, G. (2014b). Applying a resting operation to alleviate bioclogging in vertical flow constructed wetlands: An experimental lab evaluation. Journal of Environmental Management, 136, 47–53. https://doi.org/10.1016/j.jenvman.2014.01.030 Huang, X., Lei, S., Wang, G., & Zeng, B. (2020). A wetland plant, Phalaris arundinacea, accumulates nitrogen and phosphorus during senescence. Environmental Science and Pollution Research, 27(31), 38928–38936. https://doi.org/10.1007/s11356-02009285-z Huang, Y., Ortiz, L., Aguirre, P., García, J., Mujeriego, R., & Bayona, J. M. (2005). Effect of design parameters in horizontal flow constructed wetland on the behaviour of volatile 99 fatty acids and volatile alkylsulfides. Chemosphere, 59(6), 769–777. https://doi.org/10.1016/j.chemosphere.2004.11.015 Ji, B., Zhao, Y., Vymazal, J., Qiao, S., Wei, T., Li, J., & Mander, Ü. (2020). Can subsurface flow constructed wetlands be applied in cold climate regions? A review of the current knowledge. Ecological Engineering, 157, 105992. https://doi.org/10.1016/j.ecoleng.2020.105992 Jia, W., Zhang, J., Wu, J., Xie, H., & Zhang, B. (2010). Effect of intermittent operation on contaminant removal and plant growth in vertical flow constructed wetlands: A microcosm experiment. Desalination, 262(1), 202–208. https://doi.org/10.1016/j.desal.2010.06.012 Kadlec, R., & Wallace, S. (2009). Treatment Wetlands (2nd ed.). CRC Press. Karathanasis, A. D., Potter, C. L., & Coyne, M. S. (2003). Vegetation effects on fecal bacteria, BOD, and suspended solid removal in constructed wetlands treating domestic wastewater. Ecological Engineering, 20(2), 157–169. https://doi.org/10.1016/S0925-8574(03)00011-9 Karatzoglou, A., Smola, A., Hornik, K., & Zeileis, A. (2004). kernlab—An S4 Package for Kernel Methods in R. Journal of Statistical Software, 11, 1–20. https://doi.org/10.18637/jss.v011.i09 Kargol, A. K., Burrell, S. R., Chakraborty, I., & Gough, H. L. (2023). Synthetic wastewater prepared from readily available materials: Characteristics and economics. PLOS Water, 2(9), e0000178. https://doi.org/10.1371/journal.pwat.0000178 100 Kim, S.-Y., Zhou, X., Freeman, C., & Kang, H. (2022). Changing thermal sensitivity of bacterial communities and soil enzymes in a bog peat in spring, summer and autumn. Applied Soil Ecology, 173, 104382. https://doi.org/10.1016/j.apsoil.2021.104382 Knowles, P., Dotro, G., Nivala, J., & García, J. (2011). Clogging in subsurface-flow treatment wetlands: Occurrence and contributing factors. Ecological Engineering, 37(2), 99–112. https://doi.org/10.1016/j.ecoleng.2010.08.005 Kuhn, M., & Wickham, H. (2020). Tidymodels: A collection of packages for modeling and machine learning using tidyverse principles. https://www.tidymodels.org Lai, W.-L., Zhang, Y., & Chen, Z.-H. (2012). Radial oxygen loss, photosynthesis, and nutrient removal of 35 wetland plants. Ecological Engineering, 39, 24–30. https://doi.org/10.1016/j.ecoleng.2011.11.010 Lam, V. S., Tran, T. C. P., Vo, T.-D.-H., Nguyen, D. D., & Nguyen, X. C. (2024). Metaanalysis review for pilot and large-scale constructed wetlands: Design parameters, treatment performance, and influencing factors. Science of The Total Environment, 927, 172140. https://doi.org/10.1016/j.scitotenv.2024.172140 Langergraber, G. (2007). Simulation of the treatment performance of outdoor subsurface flow constructed wetlands in temperate climates. Science of The Total Environment, 380(1), 210–219. https://doi.org/10.1016/j.scitotenv.2006.10.030 Langergraber, G., Dotro, G., Nivala, J., Rizzo, A., & Stein, O. (2020). Wetland Technology: Practical Information on the Design and Application of Treatment Wetlands. IWA Publishing. https://doi.org/10.2166/9781789060171 101 Langergraber, G., Prandtstetten, C., Pressl, A., Rohrhofer, R., & Haberl, R. (2007). Removal efficiency of subsurface vertical flow constructed wetlands for different organic loads. Water Science and Technology, 56(3), 75–84. https://doi.org/10.2166/wst.2007.495 Langergraber, G., Pressl, A., Leroch, K., Rohrhofer, R., & Haberl, R. (2009). Experiences with a top layer of gravel to enhance the performance of vertical flow constructed wetlands at cold temperatures. Water Science and Technology, 59(6), 1111–1116. https://doi.org/10.2166/wst.2009.085 Langergraber, G., & Šimůnek, J. (2012). Reactive Transport Modeling of Subsurface Flow Constructed Wetlands Using the HYDRUS Wetland Module. Vadose Zone Journal, 11(2), vzj2011.0104. https://doi.org/10.2136/vzj2011.0104 Liebetrau, A. M. (1983). Measures of Association. SAGE. Lim, B., Arik, S. O., Loeff, N., & Pfister, T. (2020). Temporal Fusion Transformers for Interpretable Multi-horizon Time Series Forecasting (No. arXiv:1912.09363). arXiv. https://doi.org/10.48550/arXiv.1912.09363 Liu, H., Hu, Z., Zhang, J., Ji, M., Zhuang, L., Nie, L., & Liu, Z. (2018). Effects of solids accumulation and plant root on water flow characteristics in horizontal subsurface flow constructed wetland. Ecological Engineering, 120, 481–486. https://doi.org/10.1016/j.ecoleng.2018.07.003 Liu, S., Long, Z., Liang, J., Zhang, J., Hu, D., Hou, P., & Zhang, G. (2025). Interpretable causal machine learning optimization tool for improving efficiency of internal carbon 102 source-biological denitrification. Bioresource Technology, 416, 131787. https://doi.org/10.1016/j.biortech.2024.131787 Liu, W., Chu, Y., Tan, Q., Chen, J., Yang, L., Ma, L., Zhang, Y., Wu, Z., & He, F. (2022). Cold temperature mediated nitrate removal pathways in electrolysis-assisted constructed wetland systems under different influent C/N ratios and anode materials. Chemosphere, 295, 133867. https://doi.org/10.1016/j.chemosphere.2022.133867 Liu, Y., Han, M., Li, F., Zhang, N., Lu, S., Liu, X., & Wu, F. (2024). Performance and mechanism of SMX removal by an electrolysis-integrated ecological floating bed at low temperatures: A new perspective of plant activity, iron plaque, and microbial functions. Journal of Hazardous Materials, 463, 132802. https://doi.org/10.1016/j.jhazmat.2023.132802 López, D., Sepúlveda-Mardones, M., Ruiz-Tagle, N., Sossa, K., Uggetti, E., & Vidal, G. (2019). Potential methane production and molecular characterization of bacterial and archaeal communities in a horizontal subsurface flow constructed wetland under cold and warm seasons. Science of The Total Environment, 648, 1042–1051. https://doi.org/10.1016/j.scitotenv.2018.08.186 Luo, Y., Chen, Q., Liu, F., & Dai, C. (2023). Both species richness and growth forms affect nutrient removal in constructed wetlands: A mesocosm experiment. Frontiers in Ecology and Evolution, 11. https://www.frontiersin.org/articles/10.3389/fevo.2023.1139053 103 MacKenzie, W. H., & Moran, J. R. (2004). Wetlands of British Columbia: A guide to identification. Province of British Columbia. Martí, A. C., Pucher, B., Hernández-Crespo, C., Monerris, M. M., & Langergraber, G. (2018). Numerical simulation of vertical flow wetlands with special emphasis on treatment performance during winter. Water Science and Technology, 78(9), 2019– 2026. https://doi.org/10.2166/wst.2018.479 Mburu, N., Sanchez-Ramos, D., Rousseau, D. P. L., van Bruggen, J. J. A., Thumbi, G., Stein, O. R., Hook, P. B., & Lens, P. N. L. (2012). Simulation of carbon, nitrogen and sulphur conversion in batch-operated experimental wetland mesocosms. Ecological Engineering, 42, 304–315. https://doi.org/10.1016/j.ecoleng.2012.02.003 Mietto, A., Politeo, M., Breschigliaro, S., & Borin, M. (2015). Temperature influence on nitrogen removal in a hybrid constructed wetland system in Northern Italy. Ecological Engineering, 75, 291–302. https://doi.org/10.1016/j.ecoleng.2014.11.027 Milani, M., Marzo, A., Toscano, A., Consoli, S., Cirelli, G. L., Ventura, D., & Barbagallo, S. (2019). Evapotranspiration from Horizontal Subsurface Flow Constructed Wetlands Planted with Different Perennial Plant Species. Water, 11(10), Article 10. https://doi.org/10.3390/w11102159 Molnar, C. (2024). Interpretable Machine Learning. https://christophm.github.io/interpretable-ml-book/ 104 Molnar, C., Casalicchio, G., & Bischl, B. (2018). iml: An R package for interpretable machine learning. Journal of Open Source Software, 3(26), 786. Molnar, C., König, G., Herbinger, J., Freiesleben, T., Dandl, S., Scholbeck, C. A., Casalicchio, G., Grosse-Wentrup, M., & Bischl, B. (2022). General Pitfalls of ModelAgnostic Interpretation Methods for Machine Learning Models. In A. Holzinger, R. Goebel, R. Fong, T. Moon, K.-R. Müller, & W. Samek (Eds.), xxAI - Beyond Explainable AI: International Workshop, Held in Conjunction with ICML 2020, July 18, 2020, Vienna, Austria, Revised and Extended Papers (pp. 39–68). Springer International Publishing. https://doi.org/10.1007/978-3-031-04083-2_4 Moraffah, R., Karami, M., Guo, R., Raglin, A., & Liu, H. (2020). Causal Interpretability for Machine Learning—Problems, Methods and Evaluation. ACM SIGKDD Explorations Newsletter, 22, 18–33. https://doi.org/10.1145/3400051.3400058 Morcillo, R. J. L., & Manzanera, M. (2021). The Effects of Plant-Associated Bacterial Exopolysaccharides on Plant Abiotic Stress Tolerance. Metabolites, 11(6), Article 6. https://doi.org/10.3390/metabo11060337 Municipal Wastewater Regulation, BC Reg 87/2012 (2012). https://www.bclaws.gov.bc.ca/civix/document/id/complete/statreg/87_2012 Muñoz, P., Drizo, A., & Cully Hession, W. (2006). Flow patterns of dairy wastewater constructed wetlands in a cold climate. Water Research, 40(17), 3209–3218. https://doi.org/10.1016/j.watres.2006.06.036 105 Nedwell, D. B. (1999). Effect of low temperature on microbial growth: Lowered affinity for substrates limits growth at low temperature. FEMS Microbiology Ecology, 30(2), 101–111. https://doi.org/10.1111/j.1574-6941.1999.tb00639.x Nguyen, X. C., Nguyen, T. P., Lam, V. S., Le, P.-C., Vo, T. D. H., Hoang, T.-H. T., Chung, W. J., Chang, S. W., & Nguyen, D. D. (2024a). Estimating ammonium changes in pilot and full-scale constructed wetlands using kinetic model, linear regression, and machine learning. Science of The Total Environment, 907, 168142. https://doi.org/10.1016/j.scitotenv.2023.168142 Nguyen, X. C., Nguyen, T. P., Lam, V. S., Le, P.-C., Vo, T. D. H., Hoang, T.-H. T., Chung, W. J., Chang, S. W., & Nguyen, D. D. (2024b). Estimating ammonium changes in pilot and full-scale constructed wetlands using kinetic model, linear regression, and machine learning. Science of The Total Environment, 907, 168142. https://doi.org/10.1016/j.scitotenv.2023.168142 Nguyen, X. C., Nguyen, T. T. H., Le, Q. V., Le, P. C., Srivastav, A. L., Pham, Q. B., Nguyen, P. M., La, D. D., Rene, E. R., Ngo, H. H., Chang, S. W., & Nguyen, D. D. (2022). Developing a new approach for design support of subsurface constructed wetland using machine learning algorithms. Journal of Environmental Management, 301, 113868. https://doi.org/10.1016/j.jenvman.2021.113868 Nivala, J., Boog, J., Headley, T., Aubron, T., Wallace, S., Brix, H., Mothes, S., van Afferden, M., & Müller, R. A. (2019). Side-by-side comparison of 15 pilot-scale conventional and intensified subsurface flow wetlands for treatment of domestic wastewater. 106 Science of The Total Environment, 658, 1500–1513. https://doi.org/10.1016/j.scitotenv.2018.12.165 Nivala, J., Wallace, S., Headley, T., Kassa, K., Brix, H., van Afferden, M., & Müller, R. (2013). Oxygen transfer and consumption in subsurface flow treatment wetlands. Ecological Engineering, 61, 544–554. https://doi.org/10.1016/j.ecoleng.2012.08.028 Nivala, J., Wallace, S., van Afferden, M., & Müller, R. A. (2022). Evapotranspiration dynamics in aerated and non-aerated subsurface flow treatment wetlands. Science of The Total Environment, 843, 156605. https://doi.org/10.1016/j.scitotenv.2022.156605 Okada, K., & Vymazal, J. (2023). The effect of aboveground biomass harvesting on nutrients removal in a constructed wetland treating municipal sewage. Ecological Engineering, 190, 106918. https://doi.org/10.1016/j.ecoleng.2023.106918 Ouellet-Plamondon, C., Chazarenc, F., Comeau, Y., & Brisson, J. (2006). Artificial aeration to increase pollutant removal efficiency of constructed wetlands in cold climate. Ecological Engineering, 27(3), 258–264. https://doi.org/10.1016/j.ecoleng.2006.03.006 Paing, J., Guilbert, A., Gagnon, V., & Chazarenc, F. (2015). Effect of climate, wastewater composition, loading rates, system age and design on performances of French vertical flow constructed wetlands: A survey based on 169 full scale systems. Ecological Engineering, 80, 46–52. https://doi.org/10.1016/j.ecoleng.2014.10.029 Panasiuk, O., Hedström, A., Langeveld, J., & Viklander, M. (2022). Identifying sources of infiltration and inflow in sanitary sewers in a northern community: Comparative 107 assessment of selected methods. Water Science and Technology, 86(1), 1–16. https://doi.org/10.2166/wst.2022.151 Pinney, M. L., Westerhoff, P. K., & Baker, L. (2000). Transformations in dissolved organic carbon through constructed wetlands. Water Research, 34(6), 1897–1911. https://doi.org/10.1016/S0043-1354(99)00330-9 Põldvere, E., Karabelnik, K., Noorvee, A., Maddison, M., Nurk, K., Zaytsev, I., & Mander, Ü. (2009). Improving wastewater effluent filtration by changing flow regimes— Investigations in two cold climate pilot scale systems. Ecological Engineering, 35(2), 193–203. https://doi.org/10.1016/j.ecoleng.2008.05.019 Prost-Boucle, S., Garcia, O., & Molle, P. (2015a). French vertical-flow constructed wetlands in mountain areas: How do cold temperatures impact performances? Water Science and Technology, 71(8), 1219–1228. https://doi.org/10.2166/wst.2015.074 Prost-Boucle, S., Garcia, O., & Molle, P. (2015b). French vertical-flow constructed wetlands in mountain areas: How do cold temperatures impact performances? Water Science and Technology, 71(8), 1219–1228. https://doi.org/10.2166/wst.2015.074 Pucher, B., & Langergraber, G. (2019). Influence of design parameters on the treatment performance of VF wetlands – a simulation study. Water Science and Technology, 80(2), 265–273. https://doi.org/10.2166/wst.2019.268 Riley, K. A., Stein, O. R., & Hook, P. B. (2005). Ammonium Removal in Constructed Wetland Microcosms as Influenced by Season and Organic Carbon Load. Journal of 108 Environmental Science and Health, Part A, 40(6–7), 1109–1121. https://doi.org/10.1081/ESE-200055594 Rusten, B., Hem, L. J., & Ødegaard, H. (1995). Nitrification of municipal wastewater in moving-bed biofilm reactors. Water Environment Research, 67(1), 75–86. https://doi.org/10.2175/106143095X131213 Rusten, B., & Ødegaard, H. (2023). Nitrogen removal in moving-bed biofilm reactor plants at low temperatures: Experiences from Norway. Water Science and Technology, 87(10), 2432–2440. https://doi.org/10.2166/wst.2023.154 Salvetti, R., Azzellino, A., Canziani, R., & Bonomo, L. (2006). Effects of temperature on tertiary nitrification in moving-bed biofilm reactors. Water Research, 40(15), 2981– 2993. https://doi.org/10.1016/j.watres.2006.05.013 Samsó, R., & Garcia, J. (2013). BIO_PORE, a mathematical model to simulate biofilm growth and water quality improvement in porous media: Application and calibration for constructed wetlands. Ecological Engineering, 54, 116–127. https://doi.org/10.1016/j.ecoleng.2013.01.021 Sani, A., Scholz, M., & Bouillon, L. (2013). Seasonal assessment of experimental verticalflow constructed wetlands treating domestic wastewater. Bioresource Technology, 147, 585–596. https://doi.org/10.1016/j.biortech.2013.08.076 Sarker, I. H. (2021). Machine Learning: Algorithms, Real-World Applications and Research Directions. SN Computer Science, 2(3), 160. https://doi.org/10.1007/s42979-02100592-x 109 Seintos, T., Koukoura, A., Statiris, E., Noutsopoulos, C., Mamais, D., Μasi, F., Prado, O., Rizzo, A., Bartroli, A., Stasinakis, A. S., & Malamis, S. (2024). Long-term operation of an upflow anaerobic sludge blanket reactor coupled with a two-stage constructed wetland for domestic wastewater treatment. Chemical Engineering Journal, 157216. https://doi.org/10.1016/j.cej.2024.157216 Sharma, P. K., Takashi, I., Kato, K., Ietsugu, H., Tomita, K., & Nagasawa, T. (2013). Seasonal efficiency of a hybrid sub-surface flow constructed wetland system in treating milking parlor wastewater at northern Hokkaido. Ecological Engineering, 53, 257–266. https://doi.org/10.1016/j.ecoleng.2012.12.054 Singh, N. K., Yadav, M., Singh, V., Padhiyar, H., Kumar, V., Bhatia, S. K., & Show, P.-L. (2023). Artificial intelligence and machine learning-based monitoring and design of biological wastewater treatment systems. Bioresource Technology, 369, 128486. https://doi.org/10.1016/j.biortech.2022.128486 Singh, S., Kulshreshtha, N. M., Goyal, S., Brighu, U., Bezbaruah, A. N., & Gupta, A. B. (2022a). Performance prediction of horizontal flow constructed wetlands by employing machine learning. Journal of Water Process Engineering, 50, 103264. https://doi.org/10.1016/j.jwpe.2022.103264 Singh, S., Kulshreshtha, N. M., Goyal, S., Brighu, U., Bezbaruah, A. N., & Gupta, A. B. (2022b). Performance prediction of horizontal flow constructed wetlands by employing machine learning. Journal of Water Process Engineering, 50, 103264. https://doi.org/10.1016/j.jwpe.2022.103264 110 Soti, A., Singh, S., Verma, V., Mohan Kulshreshtha, N., Brighu, U., Kalbar, P., & Bhushan Gupta, A. (2023). Assessment of removal rate coefficient in vertical flow constructed wetland employing machine learning for low organic loaded systems. Bioresource Technology, 376, 128909. https://doi.org/10.1016/j.biortech.2023.128909 Stefanakis, A. I., & Tsihrintzis, V. A. (2012). Effects of loading, resting period, temperature, porous media, vegetation and aeration on performance of pilot-scale vertical flow constructed wetlands. Chemical Engineering Journal, 181–182, 416–430. https://doi.org/10.1016/j.cej.2011.11.108 Stein, O. R., & Hook, P. B. (2005). Temperature, Plants, and Oxygen: How Does Season Affect Constructed Wetland Performance? Journal of Environmental Science and Health, Part A, 40(6–7), 1331–1342. https://doi.org/10.1081/ESE-200055840 Stein, O. R., Hook, P. B., Biederman, J. A., Allen, W. C., & Borden, D. J. (2003). Does batch operation enhance oxidation in subsurface constructed wetlands? Water Science and Technology, 48(5), 149–156. https://doi.org/10.2166/wst.2003.0306 Sultana, M.-Y., Mourti, C., Tatoulis, T., Akratos, C. S., Tekerlekopoulou, A. G., & Vayenas, D. V. (2016). Effect of hydraulic retention time, temperature, and organic load on a horizontal subsurface flow constructed wetland treating cheese whey wastewater. Journal of Chemical Technology & Biotechnology, 91(3), 726–732. https://doi.org/10.1002/jctb.4637 111 Tanaka, N., Yutani, K., Aye, T., & Jinadasa, K. B. S. N. (2007). Effect of broken dead culms of Phragmites australis on radial oxygen loss in relation to radiation and temperature. Hydrobiologia, 583(1), 165–172. https://doi.org/10.1007/s10750-006-0483-7 Tanner, C. C., D’Eugenio, J., McBride, G. B., Sukias, J. P. S., & Thompson, K. (1999). Effect of water level fluctuation on nitrogen removal from constructed wetland mesocosms. Ecological Engineering, 12(1), 67–92. https://doi.org/10.1016/S0925-8574(98)00055X Taylor, C. R., Hook, P. B., Stein, O. R., & Zabinski, C. A. (2011). Seasonal effects of 19 plant species on COD removal in subsurface treatment wetland microcosms. Ecological Engineering, 37(5), 703–710. https://doi.org/10.1016/j.ecoleng.2010.05.007 Thorén, A.-K., Legrand, C., & Tonderski, K. S. (2004a). Temporal export of nitrogen from a constructed wetland: Influence of hydrology and senescing submerged plants. Ecological Engineering, 23(4), 233–249. https://doi.org/10.1016/j.ecoleng.2004.09.007 Thorén, A.-K., Legrand, C., & Tonderski, K. S. (2004b). Temporal export of nitrogen from a constructed wetland: Influence of hydrology and senescing submerged plants. Ecological Engineering, 23(4), 233–249. https://doi.org/10.1016/j.ecoleng.2004.09.007 Torrens, A., Molle, P., Boutin, C., & Salgot, M. (2009). Impact of design and operation variables on the performance of vertical-flow constructed wetlands and intermittent 112 sand filters treating pond effluent. Water Research, 43(7), 1851–1858. https://doi.org/10.1016/j.watres.2009.01.023 Tunçsiper, B., Drizo, A., & Twohig, E. (2015). Constructed wetlands as a potential management practice for cold climate dairy effluent treatment—VT, USA. CATENA, 135, 184–192. https://doi.org/10.1016/j.catena.2015.07.028 Ury, E. A., Arrumugam, P., Herbert, E. R., Badiou, P., Page, B., & Basu, N. B. (2023). Source or sink? Meta-analysis reveals diverging controls of phosphorus retention and release in restored and constructed wetlands. Environmental Research Letters, 18(8), 083002. https://doi.org/10.1088/1748-9326/ace6bf Varma, M., Gupta, A. K., Ghosal, P. S., & Majumder, A. (2021). A review on performance of constructed wetlands in tropical and cold climate: Insights of mechanism, role of influencing factors, and system modification in low temperature. Science of The Total Environment, 755, 142540. https://doi.org/10.1016/j.scitotenv.2020.142540 Vidal, B., Kinnunen, J., Hedström, A., Heiderscheidt, E., Rossi, P., & Herrmann, I. (2023). Treatment efficiency of package plants for on-site wastewater treatment in cold climates. Journal of Environmental Management, 342, 118214. https://doi.org/10.1016/j.jenvman.2023.118214 Vymazal, J. (2005). Horizontal sub-surface flow and hybrid constructed wetlands systems for wastewater treatment. Ecological Engineering, 25(5), 478–490. https://doi.org/10.1016/j.ecoleng.2005.07.010 113 Vymazal, J. (2007). Removal of nutrients in various types of constructed wetlands. Science of The Total Environment, 380(1), 48–65. https://doi.org/10.1016/j.scitotenv.2006.09.014 Vymazal, J. (2020). Removal of nutrients in constructed wetlands for wastewater treatment through plant harvesting – Biomass and load matter the most. Ecological Engineering, 155, 105962. https://doi.org/10.1016/j.ecoleng.2020.105962 Wallace, S., & Nivala, J. (2005). Thermal response of a horizontal subsurface flow wetland in a cold temperate climate. International Water Association’s Specialist Group on the Use of Macrophytes in Water Pollution Control, 29. Wallace, S., Parkin, G., & Cross, C. (2001). Cold climate wetlands: Design and performance. Water Science and Technology : A Journal of the International Association on Water Pollution Research, 44, 259–265. https://doi.org/10.2166/wst.2001.0838 Wang, H., Sheng, L., & Xu, J. (2021). Clogging mechanisms of constructed wetlands: A critical review. Journal of Cleaner Production, 295, 126455. https://doi.org/10.1016/j.jclepro.2021.126455 Wang, M., Zhang, D. Q., Dong, J. W., & Tan, S. K. (2017). Constructed wetlands for wastewater treatment in cold climate—A review. Journal of Environmental Sciences, 57, 293–311. https://doi.org/10.1016/j.jes.2016.12.019 Wang, Q., Xie, H., Zhang, J., Liang, S., Ngo, H. H., Guo, W., Liu, C., Zhao, C., & Li, H. (2015a). Effect of plant harvesting on the performance of constructed wetlands during 114 winter: Radial oxygen loss and microbial characteristics. Environmental Science and Pollution Research, 22(10), 7476–7484. https://doi.org/10.1007/s11356-014-3966-5 Wang, Q., Xie, H., Zhang, J., Liang, S., Ngo, H. H., Guo, W., Liu, C., Zhao, C., & Li, H. (2015b). Effect of plant harvesting on the performance of constructed wetlands during winter: Radial oxygen loss and microbial characteristics. Environmental Science and Pollution Research, 22(10), 7476–7484. https://doi.org/10.1007/s11356-014-3966-5 Wastewater Systems Effluent Regulations, SOR/2012-139 (2012). https://lawslois.justice.gc.ca/eng/regulations/sor-2012-139/fulltext.html Weerakoon, G. m. p. r., Jinadasa, K. b. s. n., Herath, G. b. b., Mowjood, M. i. m., Zhang, D., Tan, S. K., & Jern, N. W. (2016). Performance of Tropical Vertical Subsurface Flow Constructed Wetlands at Different Hydraulic Loading Rates. CLEAN – Soil, Air, Water, 44(8), 938–948. https://doi.org/10.1002/clen.201500101 Wright, M. N., & Ziegler, A. (2017). ranger: A Fast Implementation of Random Forests for High Dimensional Data in C++ and R. Journal of Statistical Software, 77, 1–17. https://doi.org/10.18637/jss.v077.i01 Xu, R., Nettleton, D., & Nordman, D. J. (2016). Case-Specific Random Forests. Journal of Computational and Graphical Statistics, 25(1), 49–65. https://doi.org/10.1080/10618600.2014.983641 Xu, X., Wang, J., Li, J., Fan, A., Zhang, Y., Xu, C., Qin, H., Mu, F., & Xu, T. (2023). Research on COD measurement method based on UV-Vis absorption spectra of 115 transmissive and reflective detection systems. Frontiers in Environmental Science, 11. https://doi.org/10.3389/fenvs.2023.1175363 Yan, Y., & Xu, J. (2014). Improving Winter Performance of Constructed Wetlands for Wastewater Treatment in Northern China: A Review. Wetlands, 34(2), 243–253. https://doi.org/10.1007/s13157-013-0444-7 Yang, B., Xiao, Z., Meng, Q., Yuan, Y., Wang, W., Wang, H., Wang, Y., & Feng, X. (2023). Deep learning-based prediction of effluent quality of a constructed wetland. Environmental Science and Ecotechnology, 13, 100207. https://doi.org/10.1016/j.ese.2022.100207 Yates, C. N., Varickanickal, J., Cousins, S., & Wootton, B. (2016). Testing the ability to enhance nitrogen removal at cold temperatures with C. aquatilis in a horizontal subsurface flow wetland system. Ecological Engineering, 94, 344–351. https://doi.org/10.1016/j.ecoleng.2016.05.064 Zhang, D. Q., Tan, S. K., Gersberg, R. M., Zhu, J., Sadreddini, S., & Li, Y. (2012). Nutrient removal in tropical subsurface flow constructed wetlands under batch and continuous flow conditions. Journal of Environmental Management, 96(1), 1–6. https://doi.org/10.1016/j.jenvman.2011.10.009 Zhang, J., Xie, H., Bai, G., Guo, D., Yang, L., Lan, J., & Ren, Y. (2024). Response characteristics of plants and pollutant removal in subsurface flow constructed wetlands under resting operation. Chemical Engineering Journal, 494, 152930. https://doi.org/10.1016/j.cej.2024.152930 116 Zhang, L., Mu, L., Xiong, Y., Xi, B., Li, G., & Li, C. (2015). The development of a natural heating technology for constructed wetlands in cold climates. Ecological Engineering, 75, 51–60. https://doi.org/10.1016/j.ecoleng.2014.11.025 Zhang, S.-Z., & Jiang, H. (2024). New Perspective for the Prediction of Pollutant Removal Efficiency in Constructed Wetlands: Using a Genetic Algorithm-Assisted Machine Learning Model. ACS ES&T Water, 4(11), 5053–5064. https://doi.org/10.1021/acsestwater.4c00635 Zhao, L., Zhu, W., & Tong, W. (2009). Clogging processes caused by biofilm growth and organic particle accumulation in lab-scale vertical flow constructed wetlands. Journal of Environmental Sciences, 21(6), 750–757. https://doi.org/10.1016/S10010742(08)62336-0 Zheng, Y., Yang, D., Dzakpasu, M., Yang, Q., Liu, Y., Zhang, H., Zhang, L., Wang, X. C., & Zhao, Y. (2020). Effects of plants competition on critical bacteria selection and pollutants dynamics in a long-term polyculture constructed wetland. Bioresource Technology, 316, 123927. https://doi.org/10.1016/j.biortech.2020.123927 Zidan, A. R. A., El-Gamal, M. M., Rashed, A. A., & El-Hady Eid, M. A. A. (2015). Wastewater treatment in horizontal subsurface flow constructed wetlands using different media (setup stage). Water Science, 29(1), 26–35. https://doi.org/10.1016/j.wsj.2015.02.003 Zuo, S., Wang, Y., Wu, J., Zhong, F., Kong, L., Chen, Y., & Cheng, S. (2024). A partial siphon operational strategy strengthens nitrogen removal performance in partially 117 saturated vertical flow constructed wetlands. Chemosphere, 142475. https://doi.org/10.1016/j.chemosphere.2024.142475 118 Appendix A: Modelling datasets Table A-1: Papers used in temperature dataset. Paper (Grebenshchykova et al., 2020) (Mietto et al., 2015) (Muñoz et al., 2006) (Nivala et al., 2019) (Sharma et al., 2013) (Wallace & Nivala, 2005) (Wallace et al., 2001) Table A-2: Papers used in treatment datasets Systems Data Min. rows Air Temp. [°C] 4 16 -10 1 37 -2.4 2 2 -8.7 12 1883 -12 3 190 -13 1 43 -17 1 25 -28 Max. Air Temp. [°C] 18 27 -8.7 27 22 21 5.3 Paper Systems Data rows HLR [L/m2/d] (Akratos & Tsihrintzis, 2007) (Allen et al., 2002) (Arias et al., 2022) (Burgos et al., 2017) (Garcı́a et al., 2003) (Grebenshchykova et al., 2020) (Hijosa-Valsero et al., 2012) (Langergraber et al., 2007) (Langergraber, 2007) (López et al., 2019) (Nivala et al., 2019) (Ouellet-Plamondon et al., 2006) (Põldvere et al., 2009) (Sani et al., 2013) (Stein et al., 2003) (Sultana et al., 2016) (Taylor et al., 2011) (Wang et al., 2015) (D. Q. Zhang et al., 2012) (Zuo et al., 2024) 5 4 1 4 8 4 2 2 1 2 12 6 1 5 3 2 20 3 6 2 15 ± 8.1 6.8 ± 0 24 ± NA 5.5 ± 0 36 ± 0 80 ± 40 50 ± 0 36 ± 11 43 ± 0 29 ± 2 74 ± 40 30 ± 0 58 ± 0.71 56 ± 17 37 ± 0 31 ± 32 3.3 ± 0 19 ± 0 37 ± 14 100 ± 0 308 5 1 12 8 16 14 9 8 4 1837 12 2 40 42 79 46 18 6 12 Min. Eff. Temp. [°C] 7 1 1 0.1 0.48 0.79 0.52 Max. Eff. Temp. [°C] 22 24 2 23 23 22 13 Effluent temperature [°C] 15 ± 5.5 8 ± 8.9 12 ± NA 13 ± 1.5 24 ± 0 15 ± 7.6 13 ± 7.7 11 ± 4.8 11 ± 4.9 16 ± 4.8 12 ± 5.2 14 ± 6.8 4.2 ± 1.6 10 ± 2.2 12 ± 6.3 20 ± 7.9 5.5 ± 4.2 5.6 ± 2.5 25 ± 0.59 9.7 ± 7.8 119 Appendix B : Supplementary model interpretation plots Figure B-1: Pearson correlation and data distribution for features in temperature dataset. X and Y axis labels are shown along the central diagonal. Figure B-2: PDP and ICE curves for the effect of loading interval [hours] on ammonia and organic matter concentrations 120 Figure B-3: PDP curves for the effect of temperature depending on the organic matter parameter for saturated systems. Figure B-4: PDP curve for the effect of temperature depending on the whether the treatment wetland is fill-and-drain (batch) operated or permanently saturated. Figure B-5: PDP curve for the effect of temperature depending on the whether the treatment wetland is fill-and-drain (batch) operated or permanently saturated. AFTER microcosm data added. 121 Figure B-6: PDP and ICE points for the effect of plants on ammonia and organic matter concentrations Figure B-7: PDP and ICE curves for the effect of recirculation rate [fraction of inflow] on ammonia and organic matter concentrations 122