# Thermal properties of treatment wetlands operated under freezing conditions

Zhanna Grebenshchykova, Nicolas Forquet, Jacques Brisson, Xavier Lachapelle-T, Yves Comeau, Florent Chazarenc

► To cite this version:

Zhanna Grebenshchykova, Nicolas Forquet, Jacques Brisson, Xavier Lachapelle-T, Yves Comeau, et al. Thermal properties of treatment wetlands operated under freezing conditions. Water Science and Technology, 2023, 88 (8), pp.2054-2067. 10.2166/wst.2023.320. hal-04564890

HAL Id: hal-04564890

https://hal.inrae.fr/hal-04564890

Submitted on 30 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

Water Science & Technology

<figure><img src="image_2.png" /></figure>

IWA
PUBLISHING

© 2023 The Authors

Water Science & Technology Vol 88 No 8, 2054 doi: 10.2166/wst.2023.320

Thermal properties of treatment wetlands operated under freezing conditions

Zhanna Grebenshchykova<sup>a,b,*</sup>, Nicolas Forquet<sup>IWA<sup>c</sup></sup>, Jacques Brisson<sup>IWA<sup>b,d</sup></sup>, Xavier Lachapelle-T<sup>e</sup>, Yves Comeau<sup>IWA<sup>f</sup></sup> and Florent Chazarenc<sup>IWA<sup>c</sup></sup>

ABSTRACT

The use of treatment wetlands (TWs) presents particular challenges in regions with sub-zero winter temperatures, due to reduced biological activity and risk of pipe breakage or clogging due to freezing. We studied the vertical temperature distribution in four pilot-scale TWs exposed to winter temperatures in order to determine the impact of operational system parameters and the role of insulation on heat conservation inside the filtering bed. The overall temperature pattern was similar in all wetlands, with a trend of increasing temperature from the surface toward the bottom during the cold season. No freezing was detected in the wetlands despite average daily temperatures as low as –20 °C. Influent water temperature and hydraulic loading had a stronger influence on TW temperatures in winter than air temperature. The vertical distribution of temperatures in TWs is more sensitive to hydraulic loading variation in the percolating operating condition than in the saturated flow with forced aeration configuration. Our results suggest that TW systems can remain operational under cold winter conditions provided the surface is properly insulated by vegetation, mulch and/or snow.

**Key words:** cold climate, insulation, thermal distribution, treatment wetland, vertical flow

HIGHLIGHTS

* In winter, temperature increases from the surface towards the bottom in vertical flow treatment wetlands.

* The vertical distribution of temperatures is more sensitive to hydraulic loading variation in the percolating operating condition than in saturated flow with forced aeration.

* Given proper surface insulation, influent water temperature and hydraulic loading have a stronger influence on treatment wetland temperatures in winter than air temperature.

* Vertical flow treatment wetlands can remain operational in winter, if the surface is properly insulated by vegetation, mulch and/or snow.

This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and
redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).

<footnote><sup>a</sup> IMT Atlantique, Bretagne-Pays de la Loire, 4 rue Alfred Kastler - La Chantrerie CS 20722, 44307 Nantes Cedex 3, France</footnote>

<footnote><sup>b</sup> Institut de recherche en biologie végétale, 4101 Sherbrooke East, Montreal, QC H1X 2B2, Canada</footnote>

<footnote><sup>c</sup> Research Unit REVERSAAL, Centre INRAE Lyon-Grenoble Auvergne-Rhône-Alpes, 5 Rue de la Doua, CS 20244, Villeurbanne Cedex 69625, France</footnote>

<footnote><sup>d</sup> Département de sciences biologiques, Université de Montréal, C.P. 6128, Succ. Centre-ville, Montreal, QC H3C 3J7, Canada</footnote>

<footnote><sup>e</sup> Ramo, 457 Rang du Ruisseau-des-Anges Sud, Saint-Roch-de-l'Achigan, QC J0K 3H0, Canada</footnote>

<footnote><sup>f</sup> Department of Civil, Geological and Mining Engineering, Polytechnique Montréal, 2900 Boulevard Edouard-Montpetit, Montréal, QC H3T 1J4, Canada</footnote>

<footnote>*Corresponding author. E-mail: zh.grebenshchykova@gmail.com</footnote>
<figure><img src="image_2.png" /></figure>

# Thermal properties of treatment wetlands operated under freezing conditions

Zhanna Grebenshchykova<sup>a,b,*</sup>, Nicolas Forquet<sup>IWA<sup>c</sup></sup>, Jacques Brisson<sup>IWA<sup>b,d</sup></sup>, Xavier Lachapelle-T<sup>e</sup>, Yves Comeau<sup>IWA<sup>f</sup></sup> and Florent Chazarenc<sup>IWA<sup>c</sup></sup>

## ABSTRACT

The use of treatment wetlands (TWs) presents particular challenges in regions with sub-zero winter temperatures, due to reduced biological activity and risk of pipe breakage or clogging due to freezing. We studied the vertical temperature distribution in four pilot-scale TWs exposed to winter temperatures in order to determine the impact of operational system parameters and the role of insulation on heat conservation inside the filtering bed. The overall temperature pattern was similar in all wetlands, with a trend of increasing temperature from the surface toward the bottom during the cold season. No freezing was detected in the wetlands despite average daily temperatures as low as –20 °C. Influent water temperature and hydraulic loading had a stronger influence on TW temperatures in winter than air temperature. The vertical distribution of temperatures in TWs is more sensitive to hydraulic loading variation in the percolating operating condition than in the saturated flow with forced aeration configuration. Our results suggest that TW systems can remain operational under cold winter conditions provided the surface is properly insulated by vegetation, mulch and/or snow.

**Key words:** cold climate, insulation, thermal distribution, treatment wetland, vertical flow

## HIGHLIGHTS

* In winter, temperature increases from the surface towards the bottom in vertical flow treatment wetlands.

* The vertical distribution of temperatures is more sensitive to hydraulic loading variation in the percolating operating condition than in saturated flow with forced aeration.

* Given proper surface insulation, influent water temperature and hydraulic loading have a stronger influence on treatment wetland temperatures in winter than air temperature.

* Vertical flow treatment wetlands can remain operational in winter, if the surface is properly insulated by vegetation, mulch and/or snow.

<footnote><sup>a</sup> IMT Atlantique, Bretagne-Pays de la Loire, 4 rue Alfred Kastler - La Chantrerie CS 20722, 44307 Nantes Cedex 3, France</footnote>

<footnote><sup>b</sup> Institut de recherche en biologie végétale, 4101 Sherbrooke East, Montreal, QC H1X 2B2, Canada</footnote>

<footnote><sup>c</sup> Research Unit REVERSAAL, Centre INRAE Lyon-Grenoble Auvergne-Rhône-Alpes, 5 Rue de la Doua, CS 20244, Villeurbanne Cedex 69625, France</footnote>

<footnote><sup>d</sup> Département de sciences biologiques, Université de Montréal, C.P. 6128, Succ. Centre-ville, Montreal, QC H3C 3J7, Canada</footnote>

<footnote><sup>e</sup> Ramo, 457 Rang du Ruisseau-des-Anges Sud, Saint-Roch-de-l'Achigan, QC J0K 3H0, Canada</footnote>

<footnote><sup>f</sup> Department of Civil, Geological and Mining Engineering, Polytechnique Montréal, 2900 Boulevard Edouard-Montpetit, Montréal, QC H3T 1J4, Canada</footnote>

<footnote>*Corresponding author. E-mail: zh.grebenshchykova@gmail.com</footnote>
## GRAPHICAL ABSTRACT

<figure><img src="image_3.png" /></figure>

# 1. INTRODUCTION

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
## GRAPHICAL ABSTRACT

<figure><img src="image_3.png" /></figure>

# 1. INTRODUCTION

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.

## 2. MATERIALS AND METHODS

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).

### 2.1. TW design

Each pilot-scale TW measured 4.5 m in length, 2.5 m in width and 1.2 m in depth. Figure 1 illustrates the filtering bed composition of the four identical pilot-scale TWs as well as the position of the temperature sensors. The beds were planted with willows (Salix miyabeana 'SX67') at a density of 4 plants/m². The system was fed with a primary effluent of municipal wastewater from the Saint-Roch-de-l'Achigan wastewater treatment facility.

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.

### 2.2. Operating conditions and treatment performance

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.

### 2.3. Temperature data collection inside the TWs

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
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.

## 2. MATERIALS AND METHODS

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).

### 2.1. TW design

Each pilot-scale TW measured 4.5 m in length, 2.5 m in width and 1.2 m in depth. Figure 1 illustrates the filtering bed composition of the four identical pilot-scale TWs as well as the position of the temperature sensors. The beds were planted with willows (Salix miyabeana 'SX67') at a density of 4 plants/m². The system was fed with a primary effluent of municipal wastewater from the Saint-Roch-de-l'Achigan wastewater treatment facility.

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.

### 2.2. Operating conditions and treatment performance

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.

### 2.3. Temperature data collection inside the TWs

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
<figure><img src="image_4.png" /><figcaption>**Figure 1** | Schematic illustration of filtering bed composition and position of the temperature sensors. TW3 and TW4 were instrumented only by sensors 0.45 and 0.75.</figcaption></figure>

<table><caption>**Table 1** | Pilot experimental conditions</caption><thead><tr><td>Pilot</td><td>Hydraulic regime</td><td>Hydraulic loading rate (m³ m<sup>-2</sup> d<sup>-1</sup>)</td></tr></thead><tbody><tr><td>TW1</td><td>Percolating</td><td>0.14</td></tr><tr><td>TW2</td><td>Saturated with aeration</td><td>0.07</td></tr><tr><td>TW3</td><td>Saturated with aeration</td><td>0.14</td></tr><tr><td>TW4</td><td>Percolating</td><td>0.07</td></tr></tbody></table>

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.

## 2.4. Meteorological data collection

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).
<figure><img src="image_4.png" /><figcaption>**Figure 1** | Schematic illustration of filtering bed composition and position of the temperature sensors. TW3 and TW4 were instrumented only by sensors 0.45 and 0.75.</figcaption></figure>

<table><caption>**Table 1** | Pilot experimental conditions</caption><thead><tr><td>Pilot</td><td>Hydraulic regime</td><td>Hydraulic loading rate (m³ m<sup>-2</sup> d<sup>-1</sup>)</td></tr></thead><tbody><tr><td>TW1</td><td>Percolating</td><td>0.14</td></tr><tr><td>TW2</td><td>Saturated with aeration</td><td>0.07</td></tr><tr><td>TW3</td><td>Saturated with aeration</td><td>0.14</td></tr><tr><td>TW4</td><td>Percolating</td><td>0.07</td></tr></tbody></table>

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.

## 2.4. Meteorological data collection

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).
## 2.5. Energy balance

The energy balance in the subsurface flow TW was determined following Kadlec & Wallace (2008) but modified to consider the input energy provided by forced aeration for TW2 and TW3 as:

$$ \text{Energy inputs} - \text{Energy outputs} = \text{Change in energy storage} \qquad (1) $$

$$ \frac{dS}{dt} = [R_N + H_a + U_{wi} + U_a] - [\lambda_m \rho E T + U_{wo} + C_{LG} + C_{LA}] \qquad (2) $$

where $dS/dt$ is the energy storage change in the wetland (MJ m⁻² d⁻¹); $R_N$ is the net radiation absorbed by the wetland (MJ m⁻² d⁻¹); $H_a$ is convective transfer from air (MJ m⁻² d⁻¹); $U_{wi}$ is energy entering with water (MJ m⁻² d⁻¹); $U_a$ is energy entering with air from forced aeration (MJ m⁻² d⁻¹); $\lambda_m$ is the latent heat of vaporization of water (MJ/kg); $\rho$ is the density of water (kg/m³); $ET$ is water lost to evapotranspiration (m/d); $U_{wo}$ is energy leaving with water (MJ m⁻² d⁻¹); $C_{LG}$ is conductive heat loss to ground (lateral and vertical) (MJ m⁻² d⁻¹); $C_{LA}$ is conductive heat loss to the air (MJ m⁻² d⁻¹) and $ET$ and $R_N$ are assumed to be negligible during the cold period due to cold air temperature and snow cover (Figure 2). Both conductive heat losses can be summed up in a single term, $C_L$, and the contribution of water can be summarized using: $U_w = U_{wi} - U_{wo}$. Therefore, Equation (2) may be rewritten as follows:

$$ -C_L = \frac{dS}{dt} - U_w - U_a \qquad (3) $$

On the right-hand side of Equation (3) are quantities that can be evaluated using the available data. The conductive heat losses ($C_L$) are then deduced by the difference. TW height can be discretized in areas of uniform properties around the sensors

<figure><img src="image_5.png" /><figcaption>**Figure 2** | Major components of the wetland energy balance in wintertime. Components in red: positive energy input, blue: energy output and gray: components negligible during wintertime (adapted from Kadlec & Wallace (2008)). A component 'energy loss with air from forced aeration' was not measured in the present study.</figcaption></figure>
## 2.5. Energy balance

The energy balance in the subsurface flow TW was determined following Kadlec & Wallace (2008) but modified to consider the input energy provided by forced aeration for TW2 and TW3 as:

$$ \text{Energy inputs} - \text{Energy outputs} = \text{Change in energy storage} \qquad (1) $$

$$ \frac{dS}{dt} = [R_N + H_a + U_{wi} + U_a] - [\lambda_m \rho E T + U_{wo} + C_{LG} + C_{LA}] \qquad (2) $$

where $dS/dt$ is the energy storage change in the wetland (MJ m⁻² d⁻¹); $R_N$ is the net radiation absorbed by the wetland (MJ m⁻² d⁻¹); $H_a$ is convective transfer from air (MJ m⁻² d⁻¹); $U_{wi}$ is energy entering with water (MJ m⁻² d⁻¹); $U_a$ is energy entering with air from forced aeration (MJ m⁻² d⁻¹); $\lambda_m$ is the latent heat of vaporization of water (MJ/kg); $\rho$ is the density of water (kg/m³); $ET$ is water lost to evapotranspiration (m/d); $U_{wo}$ is energy leaving with water (MJ m⁻² d⁻¹); $C_{LG}$ is conductive heat loss to ground (lateral and vertical) (MJ m⁻² d⁻¹); $C_{LA}$ is conductive heat loss to the air (MJ m⁻² d⁻¹) and $ET$ and $R_N$ are assumed to be negligible during the cold period due to cold air temperature and snow cover (Figure 2). Both conductive heat losses can be summed up in a single term, $C_L$, and the contribution of water can be summarized using: $U_w = U_{wi} - U_{wo}$. Therefore, Equation (2) may be rewritten as follows:

On the right-hand side of Equation (3) are quantities that can be evaluated using the available data. The conductive heat losses ($C_L$) are then deduced by the difference. TW height can be discretized in areas of uniform properties around the sensors

$$ -C_L = \frac{dS}{dt} - U_w - U_a \qquad (3) $$

<figure><img src="image_5.png" /><figcaption>**Figure 2** | Major components of the wetland energy balance in wintertime. Components in red: positive energy input, blue: energy output and gray: components negligible during wintertime (adapted from Kadlec & Wallace (2008)). A component 'energy loss with air from forced aeration' was not measured in the present study.</figcaption></figure>
and consequently the net energy change over a day (ΔS) can be expressed by:

$$ \Delta S = A \sum_{i=1}^{nb\_layer} (\rho_w c_{pw} \theta_i + \rho_a c_{pa} (\phi - \theta_i) + \rho_s c_{ps} (1 - \phi))(T_{i,j} - T_{i,j-1})h_i \quad (4) $$

where A is the TW surface (m²); ρw, ρa and ρs are densities of water, air and solids, respectively (kg/m³); cpw, cpa and cps are specific heat capacities of water, air and solids, respectively (J kg⁻¹ C⁻¹); φ is porosity (–); θi is the water content of the layer i (–); hi is the thickness of the layer i (m) and Tij is the temperature in the layer i at day j (°C).

Missing temperature values in sensor time series were imputed by linear interpolation. The net energy fluxes for the water and insufflated air can be estimated by:

$$ U_w = \rho_w c_{pw} Q(T_{in,j} - T_{out,j}) \quad (5) $$

and

$$ U_a = \rho_a c_{pa} Q_a(T_{in,j} - T_{filter,j}) \quad (6) $$

where Tin,j and Tout,j are the influent and effluent temperature at day j (°C) and Tfilter,j is the filter average temperature at day j (°C)

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

Finally, the conductive heat loss through the surface, where estimated in an attempt to distinguish it from other conductive heat losses, is determined using Equation (7):

$$ C_{LA} = A k_{g, \text{eff}} \frac{T_{\text{air}, j} - T_{\text{filter}, j}}{e_{\text{mulch}} + e_{\text{snow}}} \quad (7) $$

where kg, eff is the thermal conductivity of the insulation layer (mulch and snow) calculated using:

$$ k_{g,\text{eff}} = \frac{e_{\text{mulch}} + e_{\text{snow}}}{k_{g,\text{mulch}} + k_{g,\text{snow}}} (\text{MJ m}^{-1} \text{d}^{-1} \text{C}^{-1}) \quad (8) $$

where esnow and emulch are snow and mulch thickness, respectively (m).

## 2.6. Data, statistical analyses and modeling

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).

# 3. RESULTS

## 3.1. Meteorological conditions

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.
and consequently the net energy change over a day (ΔS) can be expressed by:

$$ \Delta S = A \sum_{i=1}^{nb\_layer} (\rho_w c_{pw} \theta_i + \rho_a c_{pa} (\phi - \theta_i) + \rho_s c_{ps} (1 - \phi)) (T_{i,j} - T_{i,j-1}) h_i \quad (4) $$

where A is the TW surface (m²); ρw, ρa and ρs are densities of water, air and solids, respectively (kg/m³); cpw, cpa and cps are specific heat capacities of water, air and solids, respectively (J kg⁻¹ C⁻¹); φ is porosity (–); θi is the water content of the layer i (–); hi is the thickness of the layer i (m) and Tij is the temperature in the layer i at day j (°C).

Missing temperature values in sensor time series were imputed by linear interpolation. The net energy fluxes for the water and insufflated air can be estimated by:

$$ U_w = \rho_w c_{pw} Q(T_{in,j} - T_{out,j}) \quad (5) $$

and

$$ U_a = \rho_a c_{pa} Q_a (T_{in,j} - T_{filter,j}) \quad (6) $$

where Tin,j and Tout,j are the influent and effluent temperature at day j (°C) and Tfilter,j is the filter average temperature at day j (°C)

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

Finally, the conductive heat loss through the surface, where estimated in an attempt to distinguish it from other conductive heat losses, is determined using Equation (7):

$$ C_{LA} = A k_{g,\text{eff}} \frac{T_{\text{air},j} - T_{\text{filter},j}}{e_{\text{mulch}} + e_{\text{snow}}} \quad (7) $$

where kg, eff is the thermal conductivity of the insulation layer (mulch and snow) calculated using:

$$ k_{g,\text{eff}} = \frac{e_{\text{mulch}} + e_{\text{snow}}}{k_{g,\text{mulch}} + k_{g,\text{snow}}} (\text{MJ m}^{-1} \text{d}^{-1} \text{C}^{-1}) \quad (8) $$

where esnow and emulch are snow and mulch thickness, respectively (m).

## 2.6. Data, statistical analyses and modeling

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).

# 3. RESULTS

## 3.1. Meteorological conditions

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.
From 12 December 2017 to 12 April 2018, the snow cover measured at the weather station was never less than 5 cm. It was above 20 cm during most of the coldest months, reaching a maximal snow cover of 55 cm in February (Figure 3). It is assumed that the TW snow cover was equal to or higher than the values recorded at the weather station (Figure 3), because the TW was more effective at trapping snow, being covered with vegetation (Kadlec 1999).

## 3.2. Influent and effluent temperatures

Average influent temperature was 8.9 °C for the cold period (from November to April) and could be considered as representative of a typical municipal influent from a short sewer network after a primary settling. At the end of December, the temperature dropped into 1.8 °C in the last tank section due to a power failure in the trailer where the settling tank was located (Figure 4).

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.

## 3.3. Filtering bed temperature in the cold period

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).

<figure><img src="image_6.png" /><figcaption>**Figure 3** | Environmental conditions during the cold period. Above: Average (red line), min and max (red ribbon) air temperature (°C). Below: Snow cover (cm).</figcaption></figure>
From 12 December 2017 to 12 April 2018, the snow cover measured at the weather station was never less than 5 cm. It was above 20 cm during most of the coldest months, reaching a maximal snow cover of 55 cm in February (Figure 3). It is assumed that the TW snow cover was equal to or higher than the values recorded at the weather station (Figure 3), because the TW was more effective at trapping snow, being covered with vegetation (Kadlec 1999).

## 3.2. Influent and effluent temperatures

Average influent temperature was 8.9 °C for the cold period (from November to April) and could be considered as representative of a typical municipal influent from a short sewer network after a primary settling. At the end of December, the temperature dropped into 1.8 °C in the last tank section due to a power failure in the trailer where the settling tank was located (Figure 4).

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.

## 3.3. Filtering bed temperature in the cold period

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).

<figure><img src="image_6.png" /><figcaption>**Figure 3** | Environmental conditions during the cold period. Above: Average (red line), min and max (red ribbon) air temperature (°C). Below: Snow cover (cm).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 4** | Evolution of mean temperature inside the filtering bed measured at depths of 0.45 and 0.75 m along with inlet/outlet temperatures, from October 2017 to April 2018.</figcaption></figure>

By comparing two depths presented in all four TWs, in the high load systems (TW1 and TW3), the difference between the coldest (0.45 m) and the warmest (0.75 m) layers was on average from 1 to 2 °C. The amplitude was larger in TWs with low hydraulic loadings (TW2 and TW4), with values ranging from 4 to 10 °C. The impact of the power failure at the end of December resulted in a net decrease of the inflow temperature and consequently of the temperatures in all beds and layers. The two saturated systems (TW2 and TW3) showed less difference in temperature between the two layers than the two percolating systems (Figure 5).

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).

## 3.4. Clustering

The temperature data at 0.45 and 0.75 m were used to create a single time series per TW (Figure 6). Air temperature and influent temperature time series were also included in the clustering analysis.

A first cut at value 6 identified the two most dissimilated groups: the TWs series on the one hand and the air temperature and influent temperature on the other. The cluster dendrogram was then cut at a vertical value of 3 to obtain four similar groups: the temperature of the two low hydraulic loading units, the temperature of the two high hydraulic loading units and the air and influent temperature (Figure 6). Inside each TW group, time series are similar with slight differences between the two types of flow conditions applied: percolating and saturated. The large decrease in the influent temperature at the end of December corresponds to the heating system outage coupled with very low outside air temperatures, together resulting in a fast drop of settling tank temperature. The temperature inside the TWs was not clearly affected by this event (Figure 7).
<figure><img src="image_7.png" /><figcaption>**Figure 4** | Evolution of mean temperature inside the filtering bed measured at depths of 0.45 and 0.75 m along with inlet/outlet temperatures, from October 2017 to April 2018.</figcaption></figure>

By comparing two depths presented in all four TWs, in the high load systems (TW1 and TW3), the difference between the coldest (0.45 m) and the warmest (0.75 m) layers was on average from 1 to 2 °C. The amplitude was larger in TWs with low hydraulic loadings (TW2 and TW4), with values ranging from 4 to 10 °C. The impact of the power failure at the end of December resulted in a net decrease of the inflow temperature and consequently of the temperatures in all beds and layers. The two saturated systems (TW2 and TW3) showed less difference in temperature between the two layers than the two percolating systems (Figure 5).

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).

## 3.4. Clustering

The temperature data at 0.45 and 0.75 m were used to create a single time series per TW (Figure 6). Air temperature and influent temperature time series were also included in the clustering analysis.

A first cut at value 6 identified the two most dissimilated groups: the TWs series on the one hand and the air temperature and influent temperature on the other. The cluster dendrogram was then cut at a vertical value of 3 to obtain four similar groups: the temperature of the two low hydraulic loading units, the temperature of the two high hydraulic loading units and the air and influent temperature (Figure 6). Inside each TW group, time series are similar with slight differences between the two types of flow conditions applied: percolating and saturated. The large decrease in the influent temperature at the end of December corresponds to the heating system outage coupled with very low outside air temperatures, together resulting in a fast drop of settling tank temperature. The temperature inside the TWs was not clearly affected by this event (Figure 7).
<figure><img src="image_8.png" /><figcaption>**Figure 5** | Temperature of the different layers inside filtering beds between November 2017 and January 2018. The 0.15 cm series is not displayed for TW1 because of the important lack of values.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6** | Resulting dendrogram after hierarchical clustering of temperature time series of the four treatment wetlands, air temperature and influent temperature. TW1: high hydraulic loading and percolating flow; TW2: low hydraulic loading and saturated flow; TW3: high hydraulic loading and saturated flow; TW4: low hydraulic loading and percolating flow.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 5** | Temperature of the different layers inside filtering beds between November 2017 and January 2018. The 0.15 cm series is not displayed for TW1 because of the important lack of values.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6** | Resulting dendrogram after hierarchical clustering of temperature time series of the four treatment wetlands, air temperature and influent temperature. TW1: high hydraulic loading and percolating flow; TW2: low hydraulic loading and saturated flow; TW3: high hydraulic loading and saturated flow; TW4: low hydraulic loading and percolating flow.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 7** | Details on four cluster series. The vertical axis represents standard normal distribution values of temperature for a cold period.</figcaption></figure>

### 3.5. Energy balance

Energy balance was calculated according to the equation presented in Section 2.5 during the cold period for all the TW systems (Table 2). The thermal losses were smaller for the percolating filter compared to saturated once fed at the same loading rate. Convective energy losses on the surface were difficult to estimate, especially regarding the impact of the insulation layer, which changed during the period considered. Energy loss through walls was also difficult to estimate, depending on the position of the TW in the system (next to another bed versus next to the outside wall insulation) (Figure 8).

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.

## 4. DISCUSSION

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.

### 4.1. Effect of air temperature

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

<table><caption>**Table 2** | Energy balance calculated in the different beds (October 2017 to April 2018)</caption><thead><tr><td>Filter</td><td>U<sub>a</sub></td><td>U<sub>w</sub></td><td>E</td><td>G</td><td>ΔS</td></tr></thead><tbody><tr><td>TW1</td><td>0</td><td>1,014</td><td>1,340</td><td>197</td><td>-328</td></tr><tr><td>TW2</td><td>60</td><td>502</td><td>931</td><td>179</td><td>-372</td></tr><tr><td>TW3</td><td>81</td><td>438</td><td>868</td><td>182</td><td>-353</td></tr><tr><td>TW4</td><td>0</td><td>615</td><td>893</td><td>150</td><td>-281</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>**Figure 7** | Details on four cluster series. The vertical axis represents standard normal distribution values of temperature for a cold period.</figcaption></figure>

### 3.5. Energy balance

Energy balance was calculated according to the equation presented in Section 2.5 during the cold period for all the TW systems (Table 2). The thermal losses were smaller for the percolating filter compared to saturated once fed at the same loading rate. Convective energy losses on the surface were difficult to estimate, especially regarding the impact of the insulation layer, which changed during the period considered. Energy loss through walls was also difficult to estimate, depending on the position of the TW in the system (next to another bed versus next to the outside wall insulation) (Figure 8).

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.

## 4. DISCUSSION

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.

### 4.1. Effect of air temperature

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

<table><caption>**Table 2** | Energy balance calculated in the different beds (October 2017 to April 2018)</caption><thead><tr><td>Filter</td><td>U<sub>a</sub></td><td>U<sub>w</sub></td><td>E</td><td>G</td><td>ΔS</td></tr></thead><tbody><tr><td>TW1</td><td>0</td><td>1,014</td><td>1,340</td><td>197</td><td>-328</td></tr><tr><td>TW2</td><td>60</td><td>502</td><td>931</td><td>179</td><td>-372</td></tr><tr><td>TW3</td><td>81</td><td>438</td><td>868</td><td>182</td><td>-353</td></tr><tr><td>TW4</td><td>0</td><td>615</td><td>893</td><td>150</td><td>-281</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>**Figure 8** | Energy balance during the study period. Red are positive values, and blue are negatives.</figcaption></figure>

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.

## 4.2. Effect of influent temperature and hydraulic loading

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 (Grebenschchykova et al. 2020).
<figure><img src="image_11.png" /><figcaption>**Figure 8** | Energy balance during the study period. Red are positive values, and blue are negatives.</figcaption></figure>

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.

## 4.2. Effect of influent temperature and hydraulic loading

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 (Grebenschchykova 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.).

### 4.3. Effect of flow operating condition

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 (Grebenschchykova 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.

### 4.4. Effect of surface insulation

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.
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.).

### 4.3. Effect of flow operating condition

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 (Grebenschchykova 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.

### 4.4. Effect of surface insulation

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.
## 5. CONCLUSION

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.

## ACKNOWLEDGEMENTS

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.

## DATA AVAILABILITY STATEMENT

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

## CONFLICT OF INTEREST

The authors declare there is no conflict.

## REFERENCES

- Cohen, J., Pfeiffer, K. & Francis, J. A. 2018 Warm Arctic episodes linked with increased frequency of extreme winter weather in the United States. *Nature Communications* **9** (1), 869.

- Environment and Climate Change Canada 2018 Historical Climate Data [WWW Document]. Available from: https://climate.weather.gc.ca/index_e.html.

- Gagnon, V., Maltais-Landry, G., Puigagut, J., Chazarenc, F. & Brisson, J. 2010 Treatment of hydroponics wastewater using constructed wetlands in winter conditions. *Water, Air, & Soil Pollution* **212**, 483–490.

- 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* **155**, 105912.

- Halekoh, U., Højsgaard, S. & Yan, J. 2006 The R package geepack for generalized estimating equations. *Journal of Statistical Software* **15** (2), 1–11.

- Hastie, T., Tibshirani, R. & Friedman, J. 2009 *The Elements of Statistical Learning: Data Mining, Inference, and Prediction*. Springer Science & Business Media, New York, NY, USA.

- Kadlec, R. H. 1999 Chemical, physical and biological cycles in treatment wetlands. *Water Science and Technology* **40** (3), 37–44.

- Kadlec, R. 2001 Thermal environments of subsurface treatment wetlands. *Water Science and Technology* **44** (11–12), 251–258.

- Kadlec, R. H. 2009 Wastewater treatment at the Houghton Lake wetland: temperatures and the energy balance. *Ecological Engineering* **35** (9), 1349–1356.

- Kadlec, R. H. & Wallace, S. 2008 *Treatment Wetlands*. CRC Press, Boca Raton, FL, USA.

- 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.

- Munoz, P., Drizo, A. & Hession, W. C. 2006 Flow patterns of dairy wastewater constructed wetlands in a cold climate. *Water Research* **40** (17), 3209–3218.
## 5. CONCLUSION

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.

## ACKNOWLEDGEMENTS

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.

## DATA AVAILABILITY STATEMENT

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

## CONFLICT OF INTEREST

The authors declare there is no conflict.

## REFERENCES

* Cohen, J., Pfeiffer, K. & Francis, J. A. 2018 Warm Arctic episodes linked with increased frequency of extreme winter weather in the United States. *Nature Communications* **9** (1), 869.

* Environment and Climate Change Canada 2018 Historical Climate Data [WWW Document]. Available from: https://climate.weather.gc.ca/index_e.html.

* Gagnon, V., Maltais-Landry, G., Puigagut, J., Chazarenc, F. & Brisson, J. 2010 Treatment of hydroponics wastewater using constructed wetlands in winter conditions. *Water, Air, & Soil Pollution* **212**, 483–490.

* 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* **155**, 105912.

* Halekoh, U., Højsgaard, S. & Yan, J. 2006 The R package geepack for generalized estimating equations. *Journal of Statistical Software* **15** (2), 1–11.

* Hastie, T., Tibshirani, R. & Friedman, J. 2009 *The Elements of Statistical Learning: Data Mining, Inference, and Prediction*. Springer Science & Business Media, New York, NY, USA.

* Kadlec, R. H. 1999 Chemical, physical and biological cycles in treatment wetlands. *Water Science and Technology* **40** (3), 37–44.

* Kadlec, R. 2001 Thermal environments of subsurface treatment wetlands. *Water Science and Technology* **44** (11–12), 251–258.

* Kadlec, R. H. 2009 Wastewater treatment at the Houghton Lake wetland: temperatures and the energy balance. *Ecological Engineering* **35** (9), 1349–1356.

* Kadlec, R. H. & Wallace, S. 2008 *Treatment Wetlands*. CRC Press, Boca Raton, FL, USA.

* 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.

* Munoz, P., Drizo, A. & Hession, W. C. 2006 Flow patterns of dairy wastewater constructed wetlands in a cold climate. *Water Research* **40** (17), 3209–3218.
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.

Prost-Boucle, S. & Molle, P. 2012 Recirculation on a single stage of vertical flow constructed wetland: treatment limits and operation modes. *Ecological Engineering* **43**, 81–84.

Prost-Boucle, S., Garcia, O. & Molle, P. 2015 French vertical-flow constructed wetlands in mountain areas: how do cold temperatures impact performances? *Water Science and Technology* **71** (8), 1219–1228.

R Core Team 2020 *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria.
Available from: https://www.R-project.org/.

Samsó, R. & García, J. 2014 The Cartridge Theory: a description of the functioning of horizontal subsurface flow constructed wetlands for wastewater treatment, based on modelling results. *Science of the Total Environment* **473**, 651–658.

Sardá-Espinosa, A. 2017 Comparing time-series clustering algorithms in r using the dtwclust package. *R Package Vignette* **12**, 41.

Smith, I. D., Bis, G. N., Lemon, E. R. & Rozema, L. R. 1997 A thermal analysis of a sub-surface, vertical flow constructed wetland. *Water Science and Technology* **35** (5), 55–62.

Wallace, S. D. & Nivala, J. A. 2005 Thermal response of a horizontal subsurface flow wetland in a cold temperate climate. *International Water Association's Specialist Group on Use of Macrophytes in Water Pollution Control* **29**, 23–30.

Wallace, S., Parkin, G. & Cross, C. 2001 Cold climate wetlands: design and performance. *Water Science and Technology* **44** (11–12), 259–265.

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.

Wittgren, H. B. & Maehlum, T. 1997 Wastewater treatment wetlands in cold climates. *Water Science and Technology* **35** (5), 45–53.

Woźniak, R., Dittmer, U. & Welker, A. 2007 Interaction of oxygen concentration and retention of pollutants in vertical flow constructed wetlands for CSO treatment. *Water Science and Technology* **56** (3), 31–38.

First received 4 July 2023; accepted in revised form 22 September 2023. Available online 4 October 2023
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.

Prost-Boucle, S. & Molle, P. 2012 Recirculation on a single stage of vertical flow constructed wetland: treatment limits and operation modes. *Ecological Engineering* **43**, 81–84.

Prost-Boucle, S., Garcia, O. & Molle, P. 2015 French vertical-flow constructed wetlands in mountain areas: how do cold temperatures impact performances? *Water Science and Technology* **71** (8), 1219–1228.

R Core Team 2020 *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing, Vienna, Austria. Available from: https://www.R-project.org/.

Samsó, R. & García, J. 2014 The Cartridge Theory: a description of the functioning of horizontal subsurface flow constructed wetlands for wastewater treatment, based on modelling results. *Science of the Total Environment* **473**, 651–658.

Sardá-Espinosa, A. 2017 Comparing time-series clustering algorithms in r using the dtwclust package. *R Package Vignette* **12**, 41.

Smith, I. D., Bis, G. N., Lemon, E. R. & Rozema, L. R. 1997 A thermal analysis of a sub-surface, vertical flow constructed wetland. *Water Science and Technology* **35** (5), 55–62.

Wallace, S. D. & Nivala, J. A. 2005 Thermal response of a horizontal subsurface flow wetland in a cold temperate climate. *International Water Association's Specialist Group on Use of Macrophytes in Water Pollution Control* **29**, 23–30.

Wallace, S., Parkin, G. & Cross, C. 2001 Cold climate wetlands: design and performance. *Water Science and Technology* **44** (11–12), 259–265.

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.

Wittgren, H. B. & Maehlum, T. 1997 Wastewater treatment wetlands in cold climates. *Water Science and Technology* **35** (5), 45–53.

Woźniak, R., Dittmer, U. & Welker, A. 2007 Interaction of oxygen concentration and retention of pollutants in vertical flow constructed wetlands for CSO treatment. *Water Science and Technology* **56** (3), 31–38.

First received 4 July 2023; accepted in revised form 22 September 2023. Available online 4 October 2023