Thermal properties of treatment wetlands operated under freezing conditions

The use of treatment wetlands (TWs) presents particular challenges in regions with sub-zero winter temperatures. Several biological and chemical processes are negatively affected by low temperatures, reducing the pollutant removal capacity of TWs (Kadlec & Wallace 2008). In addition, freezing can cause pipes to break and ice to clog the filtering bed. Subsurface flow TWs are usually preferable to other types of TWs for application in cold climate because of their better insulation (Wittgren & Maehlum 1997; Ouellet-Plamondon et al. 2006).

The majority of studies on subsurface TWs under a cold climate have focused on pollutant removal efficiency (Gagnon et al. 2010; Wang et al. 2017), with few investigating the thermal properties of the filtering bed. Wallace & Nivala (2005) studied the thermal response of a subsurface horizontal flow TW in a cold climate. They found that in summer, the highest temperature is situated closer to the surface due to important energy inputs from solar radiation, and temperature decreases with depth, while in winter, the thermal layering is inversed – with temperature increasing with depth. This important trend was observed for vertical flow (VF) TWs as well (Smith et al. 1997).

Several studies focused instead on technical adaptations of treatment design strategies (Wittgren & Maehlum 1997; Wallace et al. 2001; Prost-Boucle & Molle 2012; Prost-Boucle et al. 2015; Grebenshchykova et al. 2020). Langergraber et al. (2009) investigated the capacity of an additional top layer to improve the thermal insulation of a VF TW. The temperature of the main layer of the filtering bed with an additional 15 cm layer of gravel on top was 1–2 °C higher than in the bed without it. However, the TW with an additional top layer showed unstable pollutant removal performance due to reduced oxygen transfer inside the filtering bed. Another type of insulating layer – mulch, was studied by Wallace et al. (2001), who highlighted the positive effect of a global thermal insulation of the filtering bed, with no negative impacts on oxygen transfer (from the surface to deepest layers), plant development or pollutant removal.

The presence of aboveground biomass helps to insulate filtering beds from freezing wind events and facilitates the accumulation of snow on the TW surface. Belowground vegetation biomass may also play an important protective role under extremely cold climatic conditions, with rhizomes providing additional thermal protection against ice formation (Munoz et al. 2006).

A better understanding of factors influencing thermal distribution in filtering beds is needed in order to determine the optimal size and design configuration of VF TWs in regions under a very cold climate (Kadlec 2009).

In the present study, the vertical temperature distribution in pilot-scale TWs exposed to very low winter temperatures was monitored at different depths by temperature sensors. We aimed to (i) determine the influence of operational system parameters (hydraulic loading, flow operating condition, influent temperature) on the internal temperature of the TWs, (ii) estimate TW energy flow during the cold period and (iii) evaluate the role of surface insulation on heat conservation in the filtering bed during the colder months.

Four experimental VF pilot-scale TWs located in Saint-Roch-de-l'Achigan, Québec, Canada (45°51′29″N, 73°35′36″W) were used to evaluate different aspects of thermal vertical distribution in TWs. The climate in this region is characterized by a cold period of 5 months (from November to March) with an average temperature of −5.5 °C (extreme minimum temperature of −36.4 °C in January 2009) and a warm period from April to October with an average temperature of 14.4 °C (extreme maximum temperature of 36.1 °C in July 2018). The average annual precipitation is 1,114 mm. All average values were calculated based on data from 2009 to 2018 (Environment and Climate Change Canada 2018).

The influent was pre-treated using a primary settling tank located inside a trailer on site for the duration of the experiment and in which heaters were used, during cold periods, to maintain the temperature above 0 °C and avoid freezing of the influent. Heating was necessary to make sure that the influent temperature in the experiment is representative of that at the outlet of a settling tank dug into the ground, as they are usually set up in full-scale treatment plants.

A percolating flow mode was used for all pilot units during the summer period. During the cold period, two operating modes were tested: percolating for treatment wetlands TW1 and TW4, and saturated with artificial aeration for TW2 and TW3 (Table 1). Artificial aeration pipes were installed at a depth of 60 cm in order to oxygenate two principal gravel layers of the filtering bed without risking undesired displacement of finer sand in the deeper sand layer. The air pump was placed outside the filter zone and placed in a box (insulated and heated) to maintain a constant air temperature >3 °C during winter. Two loading rates were applied during the current experimental periods: a low load for TW2 and TW4 and a high load for TW1 and TW3 (Table 1).

The performance of the TWs was evaluated for a period of 22 consecutive months. All TWs showed high pollutant removal efficiency with a mean TSS removal of 85 ± 11% and a mean COD removal of 91 ± 4%. Ammonia nitrification was almost complete during warm periods. During cold periods, nitrification was significantly lower for TWs operated under the saturated flow condition (Grebenshchykova et al. 2020). The present study focuses on the second cold period, from October 2017 to April 2018, during which nominal loading rates were reached.

Volumetric water content and temperature sensors (5TM sensors, Decagon Devices, Inc., USA) were installed in the center of each pilot unit, at five different depths for TW1 and TW4 and at two different depths for TW2 and TW3 (Figure 1). These sensors can operate from – 40 to 50 °C. All sensors were connected to a five-channel self-contained data recorder (Decagon Em50). Recording frequency was set to 1 min. Data measured during the 2017–2018 cold season (from 27 September 2017 to 13 April 2018) are reported in this article. Temperature data for the five layers were used for the study of temperature profile presentation and energy balance calculation. The influent temperature was measured at least twice per month (samples analyzed by a temperature laboratory sensor on-site) in the last compartment of the primary settling tank from which the wastewater was drawn to supply the TWs.

Average, minimum and maximum daily air temperature during the study were provided by the L'Assomption meteorological station, located 13 km from the experimental site (Environment Canada station #7014160; 45°48′34″N, 73°26′05″O, 21 m). The experimental site and the meteorological stations of L'Assomption and Saint-Roch-de-L'Achigan are all located in a flat plain, which minimizes micro- and mezoclimatic differences and allows the assumption that temperature, humidity and snow cover observations are similar in the three sites. The snow cover data was provided by the Saint-Roch-de-l'Achigan station (Environment Canada station #7017698; 45°50′57″N, 73°41′29″O, 57 m).

Equation (6) represents the potential heat transfer and not the real one. The incoming air temperature is considered equal to the influent temperature.

The daily temperatures inside the different layers were compared between the different beds' conditions and air temperature during the study period. Hierarchical clustering was chosen as a way to compare temperature time series of all TWs by computing the average temperature at every time step (Sardá-Espinosa 2017). The dtwclust package was used to produce a dendrogram – a binary tree where the height of each node is proportional to the value of the inter-group dissimilarity between its two daughter nodes (Halekoh et al. 2006; Hastie et al. 2009). The chosen metric for distance is the Dynamic Time Wrap (DTW) – specifically developed for time series. It is a dynamic programming algorithm that compares two series and tries to find the optimum warping path between them under certain constraints (Sardá-Espinosa 2017). Once the clusters were obtained, different cluster groups were visualized to focus on dissimilarities between them. Energy flows were estimated for a cold period according to Figure 2. The energy lost by water was estimated for all four bed conditions. The statistical analyses were carried out using R (R Core Team 2020).

The cold season, from 29 September 2017 to 13 April 2018, was characterized by an average air temperature of −1.4 °C, a minimum of −32.5 °C and a maximum of 18.8 °C.

The difference between the influent and effluent temperatures during the cold period was always less than 2 °C (Figure 4). Therefore, the system never clogged as a result of freezing.

First, we focus on the temperature recorded inside all the filters at depths of 0.45 and 0.75 m, along with influent and effluent temperature, to highlight a clear tendency in the different filters without providing too much visual information (Figure 4).

The overall temperature pattern was similar in all beds, with a trend of temperature increasing from the surface toward the bottom during the cold season. At the beginning of November, the outside air temperature dropped rapidly below 0 °C and the snow cover was not yet present on the surface (Figure 3). This resulted in the temperatures inside TWs dropping below 10 °C. From December to March, the outside air temperature dropped and remained around −10 to −20 °C. During this period, the snow cover was permanently present at the surface of the filters, with depths ranging from 20 to 50 cm. The temperature pattern was well-balanced for all filters, with the inside temperature following the same trend as the influent temperature (Figure 4).

Temperatures were always warmer and more uniform in the percolating system fed with the higher load (TW1), whereas in the three other beds, the vertical gradient (coldest temperature in surface and highest in deepest layer) was more pronounced. The largest amplitude between different layers was recorded in the percolating filter with the lowest hydraulic loading rate (TW4), with the temperature in bottom layers below 5 °C during most of the cold period (Figure 5).

In October 2017, a net energy deficit occurred in all beds, which corresponded to a temperature decrease and the absence of an insulating layer of snow. A larger energy deficit was calculated during most of the cold period, the amplitude of energy benefits and deficits being greater in highly loaded TW1 and TW3 as a result of the higher energy input from the influent.

The temperature of the different layers inside the TWs and the effluent temperature in winter were mostly influenced by the air temperature, snow cover, operating conditions (influent temperature, hydraulic (organic) loadings and flow mode) and system design. The discussion is thus structured according to each of these parameters.

Sub-zero air temperature for a period of 3–5 months per year represents the main challenge to ensuring successful operation of TW systems in a cold climate. During the present study, winter air temperatures dropped to nearly −30 °C. Such cold winter temperatures can alternate with relatively warm periods, during which mean daily temperatures can rise above 0 °C, leading to snowmelt and ice formation, altering the insulating layer formed by the snow cover and influencing the temperature inside the filtering bed. During the cold period, daily air temperature rose above 0 °C several times (Figure 3). Such great fluctuations in winter temperature are expected to increase in the future with global climate change (Cohen et al. 2018). In this study, surface insulation conditions were sufficient to avoid snowmelt events.

The high thermal inertia inside the bed linked to the presence of a snow cover has rarely been reported in the literature. During our study period, the presence of a snow cover prevented low air temperatures from having a significant impact on the filtering bed temperature. Implementing design parameters that favor the presence and establishment of a deep snow cover can therefore be highly beneficial.

The energy balance shows that the influent temperature is the main source of thermal energy for the TW during the cold period (Figure 2). Yet, influent temperature has rarely been the focus, or even considered, in TW research studies to date. While it may be an external parameter because it is generally not an option in full-scale applications, our results and other available data suggest that an influent temperature above 5 °C prevents freezing even under harsh winter conditions. Thus, the cold influent load carried by sewer networks that feed the wastewater treatment plant could have a deleterious impact on the functionality of a TW system.

Hydraulic loading is the parameter with the second greatest influence on TW temperature. As shown by the hierarchical clustering analyses, the two different levels of hydraulic loadings applied during the cold period – 0.07 and 0.14 m³ m⁻²d⁻¹, had a greater impact on bed temperature than the different flow operating conditions.

However, the estimated energy gains for the two TWs with different hydraulic loadings showed two almost identical responses. The probable hypothesis to explain this absence of difference is the low range of hydraulic loadings. The treatment efficiency of all TWs tested was high, with the few differences noted mostly due to the flow operating conditions rather than the difference in hydraulic loadings applied (Grebenshchykova et al. 2020).

Generally, the strategy adopted to maintain a positive winter energy balance in TWs is to prioritize actions that will minimize energy loss flows (Wallace & Nivala 2005). This involves estimating and controlling energy gains and losses to maintain a ‘balance point’ between energy flows. Although a potential increase in hydraulic loadings will not significantly modify the capacity of a TW to store energy, such an action may be necessary to reduce temperature differences between different deeper layers. Ensuring a more stable temperature throughout the filtering bed allows diverse microorganism pollutant removal processes to continue even if a gradient a gradient in microorganism quantity and diversity will remain due to higher pollutant concentration near the inlet (Woźniak et al. 2007; Samsó & García 2014). The objective guiding such a strategy is to determine the best compromise between, on the one hand, applying the largest organic load possible in order to enable satisfactory treatment performance, and, on the other hand, the largest hydraulic loading rate possible, to maximize the amount of energy in the wetland without generating hydraulic failure (overload, clogging, etc.).

It is generally assumed that subsurface saturated flow is the best design for cold climates (Wang et al. 2017), reducing contact with the atmosphere as well as limiting energy transfer between gas and liquid. This design parameter has an only slight offset between the two temperature series in our study. The flow operating conditions had no significant influence on TW temperature. This is most likely a consequence of the low hydraulic loadings applied during the experiment (Grebenshchykova et al. 2020).

A large batch loading generally pushes the air that is inside the wetland (warm) through the ventilation pipes, replacing it with (very cold) outside air. At a low hydraulic loading rate and therefore with few batch events, it can be hypothesized that warm air remains inside the wetland system, and that the water percolates and exchanges its energy with minimal outside air.

The hottest temperature inside the sand layer for saturated TWs could be a main result of the depth position and could be a result of smaller grain size for percolating TWs.

The natural and low-cost materials used for surface insulation in this study – a willow mulch layer of 15 cm and vegetation (Salix miyabeana), were sufficient to trap enough snow to protect the filtering beds from freezing during the cold period. The high thermal inertia of the TW units was demonstrated by a period of 7 days without feeding (without influent energy input) when air temperatures were extremely low (below −20 °C). Under these conditions, the measured TW temperatures remained above 5 °C. This response shows important potential for full-scale application.

In the present study, tall stems (more than 1.5 m) of Salix miyabeana planted at high density (4 plants/m²) acted a very effective snow trap (see Supplementary Information).

Vegetation-based insulation must be coupled with another type of insulator during the initial years of operation or when the species used provides too little biomass to protect the system from early freezing (Wittgren & Maehlum 1997). Frequently used in cold climatic conditions and recommended by previous research (Wallace et al. 2001), a mulch layer of willow wood chips was the principal insulation layer used in the present study. Mulch quality can strongly influence its thermal efficiency, and the willow wood chips have been shown to play an important role in insulation efficiency (10 cm of mulch has a thermal conductivity of 0.0052 MJ m⁻¹ d⁻¹ °C⁻¹), providing a ‘safety layer’ when snow cover is lacking, especially at the beginning of the cold period. This value is two times smaller than the thermal conductivity of the snow (0.01 MJ m⁻¹ d⁻¹ °C⁻¹) and five times smaller than the one of dry gravel (0.026 MJ m⁻¹ d⁻¹ °C⁻¹) (Kadlec 2001; Kadlec & Wallace 2008).

Indeed, the snow cover is a significant insulator of TWs (Kadlec 2001; Prost-Boucle et al. 2015). Although Wittgren & Maehlum (1997) suggest that snow cover cannot be taken for granted as an insulating layer because of its direct dependency on weather conditions, TWs must be designed to maximize surface snow. In the present study, the conditions for fast and durable snow accumulation (snow cover above 20 cm for a period of three consecutive months) were met. Among the most important factors ensuring good snow insulation are the presence of vegetation, installation of distribution pipes under the mulch layer and a saturation level below the surface. Creating a ‘dry’ layer of mulch, as typically recommended for horizontal flow TWs, and avoiding open water on the surface help build effective long-term insulation (Wallace et al. 2001; Kadlec & Wallace 2008).

In early November, when the daily air temperature dropped below 0 °C, the snow cover began to be continuously present on the TW units. From then on, the temperature inside the TWs stopped falling and remained stable throughout the rest of the cold period. This confirms the significant role played by the snow layer in insulating TWs.

Lower air temperature resulted in a decrease in water temperature inside the filtering bed at the beginning of the cold season, but surface insulation protected the beds from further temperature decreases later in winter. Two external parameters, influent temperature and hydraulic loading, became more important in influencing TW temperature in winter. The first is rarely considered in the configuration and operation of TWs. More studies are needed to fully understand the effects of this parameter.

Proper surface insulation (including vegetation, mulch and snow) can prevent hydraulic failure due to freezing of the filtering bed. Efficient insulation can generate high thermal inertia, ensuring that TW systems continue to operate under cold conditions without additional technical components.

Higher hydraulic loadings resulted in more homogeneous temperatures between the TW layers. The vertical distribution of temperatures in TWs is more sensitive to hydraulic loading variation in the percolating operating condition than under saturated flow with forced aeration.

Future studies are needed to understand the horizontal temperature distribution in full-scale TWs.

We are grateful to all of the partners who contributed financially to the PhytoValP project: Ramo, Bionest, Arcelor Mittal, Minéraux Harsco, Équipe Indigo, Ressources Aquatiques Québec, the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Consortium de recherche et d'innovation en bioprocédés industriels du Québec (CRIBIQ). Our sincere thanks also to the following partners: Naturally Wallace Consulting, SINT, GHD, NordikEau, SÉPAQ, INRAE (formerly IRSTEA), MDDLCC, MAPAQ and the Municipality of Saint-Roch-de-l'Achigan for access to its WRRF. We are also grateful to Pascale Mazerolle and Roselyne Gagné-Turcotte for field assistance and to Karen Grislis for linguistic revision of the manuscript.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.