Plant-wide modelling and analysis of WWTP temperature dynamics for sustainable heat recovery from wastewater


 Wastewater heat recovery upstream of wastewater treatment plants (WWTPs) poses a risk to treatment performance, i.e. the biological processes. In order to perform a sustainability analysis, a detailed prediction of the temperature dynamics over the WWTP is needed. A comprehensive set of heat balance equations was included in a plant-wide process model and validated for the WWTP in Linköping, Sweden, to predict temperature variations over the whole year in a temperate climate. A detailed model for the excess heat generation of biological processes was developed. The annual average temperature change from influent to effluent was 0.78 °C with clear seasonal variations; 45% of the temperature change arises from processes other than the activated sludge unit. Hence, plant-wide energy modelling was necessary to predict in-tank temperature in the biological treatment steps. The energy processes with the largest energy gains were solar radiation and biological processes, while the largest losses were from conduction, convection and atmospheric radiation. Tanks with large surface areas have a significant impact on the heat balance regardless of biological processes. Simulating a 3 °C lower influent temperature, the temperature in the activated sludge unit dropped by 2.8 °C, which had a negative impact on nitrogen removal.


INTRODUCTION
As resource recovery from wastewater (WW) is widely applied to improve sustainability of wastewater treatment (WWT) systems (Verstraete et al. 2009), wastewater heat recovery (WWHR) is gaining wider interest (Culha et al. 2015;Cecconet et al. 2020;Wärff et al. 2020). The energy in the WW largely comes from domestic hot water. According to several international studies, the energy use in the urban water cycle adds up to 10% or more of the total national energy use (Olsson 2012). Out of this only about 10% (corresponding to 1% of the total energy use) is used for withdrawal, treatment and distribution of tap water and collection, treatment and discharge of wastewater. The remaining partthe absolute majorityis used by the customers mainly for heating tap water for showers, dishwashers and laundry (Olsson 2012). Swedish Energy Agency (2009) has estimated the heat requirement in households for domestic hot water (DHW) to be 1,150 kWh/cap/yr. This contributes to an elevated temperature of wastewater compared to tap water. Wastewater heat recovery (WWHR) is a known practice in many countries. It can be implemented either up-stream of the wastewater treatment plant (WWTP), such as in showers, in building sewer stems, in sewer mains, or at the WWTP effluent after treatment (Arnell et al. 2017). From an energy recovery point of view, the more up-stream locations are more favourable (i.e. they have a higher WW temperature). On the other hand, more downstream locations have larger and more stable flows while the temperature at the same time is lower due to energy loss in the sewer network and mixing with intrusive water of lower temperature. This, along with maintenance concerns, makes passive heat exchangers more common in appliances and buildings while heat pumps are used in precincts or at WWTP effluents (Culha et al. 2015;Cecconet et al. 2020;Jonsson et al. 2020). Recovering heat from WW up-stream of the WWTP will impact the influent temperature of the WW. Preliminary studies indicate that this effect is in the range of 1-3°C (Arnell et al. 2017).
Temperature is well known to impact performance and efficiency of WWT processes. Most processes are somewhat temperature dependant with density, viscosity and other basic characteristics of fluids varying with temperature. However, biological processes are the most sensitive within the temperature range WWTPs operate in temperate climates, 6-27°C (la Cour Jansen et al. 1992;Henze et al. 2002). Rates of biological processes are reduced as temperature drops. At WWTPs with biological nitrogen removal, nitrification is the rate limiting process (Henze et al. 2008). From 20 to 10°C, the maximum (i.e. non-limited by for example dissolved oxygen concentration) nitrification rate is reduced by more than 50% (Henze et al. 2002;Gernaey et al. 2014). Consequently, there is a risk that up-stream recovery of heat from wastewater deteriorates the nitrogen removal at the WWTP and leads to higher discharges of nitrogen compounds to the environment or a higher consumption of energy and other resources at the plant to maintain the same effluent quality. The risk will depend on local conditions, such as size and quality of the sewer network and WW load relative to the design capacity of the plant. Therefore, the risk must be evaluated case by case before allowing WWHR in the catchment. To assess the plant performance as a function of temperature, process modelling and simulation can be used. Plant-wide modelling has been shown to be instructive for analysing energy use and recovery at WWTPs (Fernandez-Arevalo et al. 2017;Zaborowska et al. 2021). The bioprocess models used in research and industry, e.g. the Activated Sludge Models (ASM1, 2, 3) (Henze et al. 2000) or in commercial simulation platforms, include temperature dependency for the biological reactions.
However, the temperature across the WWTP is not constant. Solar radiation, wind, ambient temperature etc. can heat or cool uncovered tanks at WWTPs and studies have shown a temperature change of 1-2°C from influent to effluent (la Cour Jansen et al. 1992;Makinia et al. 2005;Fernandez-Arevalo et al. 2014). In literature, several examples of how to model temperature, i.e. perform an energy balance, in WWTP tanks have been presented. In a pioneering work, la Cour Jansen et al. (1992) assessed the problem of operating nitrifying WWTPs in a cold climate. A temperature model for an activated sludge tank was presented based on a simple energy balance including the significant energy influences. The model was integrated with an ASM1 process model and validated against a one-year data set from one Danish WWTP. Sedory & Stenstrom (1995) presented a temperature model for an aeration basin and validated it for five North American WWTPs. Thirty-day temperature series for different temperature intervals were used as inputs to the model. Temperatures ranged from 12 to 31°C and their model did not consider effects on WWT performance. Their model accuracy was in the range 2-3°C and did not predict extreme variations very well. In a more recent paper, Makinia et al. (2005) developed a temperature model for an activated sludge tank and compared a partial differential equation formulation with the more common 'tanks in series' approach. The model was tested for a short (a few days) data series with a temperature around 19°C and did not consider WWTP performance. Lippi et al. (2009) sought to improve the temperature model from Sedory & Stenstrom (1995). The heat flux equations were updated to more detailed alternatives. When the model was tested against the same five data sets and a few more, the model accuracy was improved. However, the model still only covered an aeration basin and did not consider the impact on other treatment processes.
Considering two facts: (1) heat recovery is of most interest in areas with a cold climate and (2) the impact of WWHR on influent temperature is in the same range as the temperature change over the WWTP, a plant-wide process model including heat balances for temperature prediction is needed. The model must integrate process performance and temperature prediction and be proven valid for relevant climate and temperature conditions. Simulations for extended (full year) time periods are required to assess the impact and potential of WWHR in WW systems. Such a model has previously not been presented in literature.
For this study, a novel plant-wide dynamic temperature model has been developed and implemented in a general process model framework. Modelled processes include grit chambers, pre-aeration tanks, primary clarifiers, activated sludge tanks, secondary settlers, tertiary biological and polishing (final clarifiers) steps and digesters. The temperature model has been developed and validated for a full year for a location with large temperature and climate variations in a temperate climate zone (Linköping, Sweden). The model was used to assess the impact of temperature change due to up-stream wastewater heat recovery in various scenarios.

General plant-wide temperature model
To model temperature over the plant and in each sub-process, Equation (1) was used. The net heat flux (f n ) must be modelled including all relevant energy gains and losses. The heat flux equation (Equation (2)) was used as a basis for the model (Makinia et al. 2005), The components of the heat flux are the following: f sr is solar radiation; f ar is atmospheric radiation; f c is conduction and convection; f e is evaporation; f a ¼ f as þ f al is aeration, sensible and latent heat losses; f m is mechanical energy (not included); and f bp represents biological processes.
Mathematical relations describing the included heat flux components were developed and selected based on the model objective, i.e. plant-wide temperature modelling for full year simulations of a WWTP in a temperate climate. The large span in water and air temperatures forced the use of dynamic expressions for most parameters. The full list of parameters is presented in Table 1  The short-wave net radiation (solar radiation, f sr [J/d]) can be calculated from direct measurements of solar radiation using a local weather station or collected from a meteorological measuring station sufficiently close to the plant as: where, SR is the solar radiation [W/m 2 ]; A s is the tank area exposed to the atmosphere [m 2 ]; and 86,400 is a time conversion factor [s/d].
The long wave net radiation (atmospheric radiation, f ar [J/d]) is based on the Stefan-Boltzmann law and is given by the following expression, where, 1 ar is the emissivity of the water surface [-]; s is the Stefan Boltzmann constant [W/(m 2 K 4 )]; l ar is the water surface reflectivity [-]; T Ã a is the absolute temperature of the ambient air [K]; T Ã w is the absolute temperature of the water surface [K]; and b ar is the atmospheric radiation factor [-].
Surface convection (f c [J/d]) can be modelled as a function of wind speed and the temperature difference between the mass of water and the air above it according to Lippi et al. (2009): where, h is the heat transfer coefficient [W/(m 2 K)] and T Ã w is the absolute temperature of the water surface [K]. For the presented model, a temperature dependent parameter expression has been adopted and implemented for Equation (5) in addition to the description by Lippi et al. (2009), see Table 1.
For the heat flux for evaporation (f e [J/d]) Equation (6) was used (Lippi et al. 2009): L c is dependent on tank size and geometry. Nu Dynamic expressions for m vap adopted from (Huber et al. 2009). See S.I. for full expressions and parameters for m vap and r vap .
Regression model adopted from Nellis & Klein (2012), which is fitted to data by Bolz & Tuve (1976).  Henze et al. 2002) See where, h m is the mass transfer coefficient [m/s]; h fg is the latent heat of vaporization of water [J/kg]; r v,Ta is the saturated vapour density at the ambient air temperature [kg/m 3 ]; and r v,Tw is the saturated vapour density at the water surface temperature [kg/m 3 ]. The applicability of Equation (6) was improved by adopting dynamic expressions for essential parameters, see Table 1.
Heat loss from aeration (f a [J/d]) is modelled as the sum of sensible heat loss (f as [J/d]) and latent heat loss (f al [J/d]), which is given by the following expressions (Makinia et al. 2005): The biological heat production from the degradation of organic matter (oxidation of carbon), nitrification and denitrification is calculated from the enthalpy of reaction (DH) of the chemical reactions (Henze et al. 2008;Kleerebezem & van Loosdrecht 2010). Dynamically calculated conversion rates in the process model were used.
where, D COD aer is the amount of COD that is used for aerobic respiration per day [g COD/d]; DN nit is the amount of nitrogen that is nitrified per day [gN/d]; DN denit is the amount of nitrogen that is denitrified per day [gN/d]; DH aer is the enthalpy of reaction for aerobic respiration [J/gCOD]; DH nit is the enthalpy of reaction for nitrification [J/gN]; DH denit is the enthalpy of reaction for denitrification [J/gN]; b aer is the metabolic efficiency for COD respiration [-]; b nit is the metabolic efficiency for nitrification [-]; b denit is the metabolic efficiency for denitrification [-]. The enthalpy of reaction for aerobic respiration is not trivial to calculate theoretically and was adopted from Blackburn & Cheng (2005). The enthalpies of reaction for nitrification and denitrification can be calculated by combining the half reactions for each conversion step using tabulated numbers for their respective enthalpy of formation (Kleerebezem & van Loosdrecht 2010;Fernandez-Arevalo et al. 2014), see Supplementary Information for reactions and calculation of enthalpies. The included enthalpies describe the total free energy from the reactions of the substrates, i.e. both the enthalpy bound in biomass and enthalpy lost as heat (Henze et al. 2002(Henze et al. , 2008Griffiths 2012). The metabolic efficiency factors (b) for each of the processes was used to compensate (Henze et al. 2002).

Plant-wide model integration
For general purpose use, the developed temperature model was implemented in the Benchmark Simulation Model No. 2 (BSM2). The BSM2 is a modelling framework for benchmarking of process control using a standardised model set-up, plant layout, load profiles and simulation and evaluation procedure (Gernaey et al. 2014). The evaluation criteria include major cost items (e.g. aeration energy and carbon source addition) in an operational cost index.
The model code for the temperature model in BSM2 is made available through Github. See S.I. for further information.

Case study Linköping wastewater treatment plant
The Linköping WWTP treats municipal and industrial wastewater (approximately 20% of annual organic load is industrial) from a total of 180,000 PE. The treatment train has primary treatment with grit removal, pre-aeration and chemically enhanced primary clarifiers; secondary biological treatment in three parallel activated sludge units (ASUs) with intermittent aeration; and finally, an advanced tertiary treatment with ozone for micro-pollutant degradation, post nitrification and denitrification in a moving bed biofilm reactor (MBBR) and post-precipitation and clarification. Mixed sludge is stabilized through a belt thickener and anaerobic digestion in three, interconnected digesters (i.e. with circular mixing inbetween). The sludge is dewatered in screw presses and the supernatant is treated in a side stream SHARON process for nitrogen removal (Figure 1). All main tanks are uncovered. A plant-wide process model for the WWTP was used for implementation and validation of the temperature model. The model framework from BSM2 was used as a basis for the model development (Gernaey et al. 2014). For the model set-up, the primary settler is modelled using an empirical mass balance model by Otterpohl & Freund (1992). A modified version of ASM1 was chosen for the activated sludge process (Gernaey et al. 2014) followed by a point-settler (with an added sludge blanket reactor volume to include settler denitrification) model for the secondary clarifier. For the ozonation, only COD degradation and oxygen addition was modelled, no micro-pollutant removal. The post nitrification process was modelled with a simplistic MBBR-model assuming a fixed surface-based nitrification rate on the carriers limited by substrate and dissolved oxygen (DO) availability and adjusted for temperature. Post denitrification was modelled using the concept of apparent kinetics (Baeten et al. 2019), using a modified ASM1 where heterotrophic microorganisms grow onand consequently are limited byprimarily NH 4 -N and secondly NO 3 -N using a switching function (Hiatt & Grady 2008). The clarifiers in the tertiary treatment were only included for temperature modelling. In the sludge line, the digestion process is modelled using ADM1 and the thickener and dewatering units are described by simple ideal separation models (Gernaey et al. 2014). The AD is heated with district heating to a constant 37°C. The side stream nitrogen removal was simplistically modelled using fixed effluent levels and removal factors for the different components based on empirical relations with data. The model was rigorously calibrated and validated following established principles for characterization (Melcer et al. 2003) and the Good Modelling Practice scheme (Rieger et al. 2012) using two separate data sets. The model was implemented, and simulations were run in Matlab/Simulink (Matlab 2017b, The Mathworks Inc., Natwick, MA, USA, 2017).
At the plant, extensive data collection was first carried out for input characterization, calibration and validation of both the regular process model and the temperature model. Along with information about the plant (tank dimensions, machine equipment, control strategies, etc.) analytical and operational data was collected for calibration (1 April-1 May 2019) and validation (19 Sept. 2018-19 Sept. 2019. Secondly, for the temperature model, meteorological data (air temperature, relative humidity, wind speed, air pressure, rain and solar radiation) was recorded for the period 30 April 2019-28 April 2020 using a Davis Instruments Vantage Pro2™ (#6162-EU) wireless weather station installed at the wastewater treatment plant in Linköping. Temperature was measured for the whole or part of the period at the influent, ASU tanks, secondary effluent (part of time), ozonation effluent (part of time) and plant effluent. The measurement uncertainty for all water temperature measurements except for the influent was +0.2°C. The influent temperature measurement uncertainty is not assessed but according to manufacturer specifications it should not exceed +0.6°C. The plant effluent temperature data had a period of faulty measurements when the sensor was periodically out of the water (23 Aug.-26 Aug. 2019). This data was excluded from the evaluation. Meteorological data was lost for three shorter periods (9 d in June/July 2019, 2 d in Feb. and 1 h in March 2020) due to limited storage capacity in the logger. For those periods data from the most near-by sites of the Swedish Meteorological and Hydrological Institute (SMHI) was used to fill the gaps: station Malmslätt, Linköping for all meteorological parameters except solar radiation (approximately. 10 km from WWTP) and station SMHI headquarters, Norrköping for solar radiation (approximately. 40 km from WWTP). Uncorrected Proof The process model was calibrated following the procedures and stop criteria in the Good Modelling Practice protocol (Rieger et al. 2012). No calibration was performed for the temperature model in the sense of adjusting model parameters to improve the fit of the model. The resulting temperature profile for effluent water temperature of a full year simulation was compared to data for validation of the model. Goodness of fit of the simulated effluent temperature to data for the validation period was calculated using Mean Absolute Error (MAE), Root Mean Squared Error (RMSE) and Square of Residuals (R 2 ), as presented in Equations (11)-(13): where, y is the simulated data; d is the measured data; N is the number of observations; and d is the average of the measured data.
Considering the intended use of the modelto evaluate the impact of influent temperature changes in the range 1-3°Cthe target for model accuracy (MAE) was set to 0.5°C in the plant effluent.

Simulation scenarios
The impact of WWHR up-stream the WWTP was simulated by imposing a reduced influent temperature with 1, 2 and 3°C in a series of three simulations (Scenarios 1, 2 and 3) to be compared with the validated reference case (Default case). Primarily, the nitrogen removal was evaluated in the activated sludge system and nitrogen concentrations at the plant effluent were compared to the effluent standards permitted.

Plant-Wide analysis of heat balance
The most relevant ambient meteorological data is shown in Figure 2. The summer period was relatively normal but the winter 2019/2020 was unusually warm with only very limited periods of really cold weather, i.e. with air temperatures below 0°C (Figure 2(a)). As the location is relatively far to the north, the solar radiation is low or negligible during winter (Figure 2(b)). For this study, the meteorological data was collected on-site, providing the most accurate data for simulations. For example, solar radiation has a high local variability depending on cloud cover. Previous studies mostly used data from off-site meteorological stations, often many kilometres away from the plant (Sedory & Stenstrom 1995;Makinia et al. 2005;Lippi et al. 2009).
The overall prediction of the temperature change over the WWTP for the full-year simulation is plotted in Figure 3(a). Comparing the simulated effluent temperature with the measured values showed a good fit. The mean absolute error was 0.30°C, which is below the target of + 0.5°C. The RMSE equals 0.41°C and the overall coefficient of determination R 2 0.98. Minor deviations can be observed in Figure 3. There was a slight over-prediction during the summer (Figure 3(b)), while the fit is generally better for the winter period (Figure 3(c)). At the end of the simulated period, as the temperature rises, model over-prediction can be seen again. During model set-up, some fixed parameters used in previous studies (Makinia et al. 2005;Lippi et al. 2009) were exchanged for dynamic expressions (e.g. r air , m air and D AB ) to provide sufficient fit over the whole temperature range. This proves the necessity to validate the model for the full range of climate conditions for which it will be used to assure predictive power. A comparison of the modelled temperature in the ASU tanks with measured values in the same position is shown in Figure 4. The fit is good also in the ASU, however, the data show faulty behaviour in Figure 4(b) with frozen values followed by negative spikes at certain periods. The measurement uncertainty of +0.2°C for the temperature sensors down-stream the influent was smaller than the model error and thus sufficiently accurate for the model validation. The relatively larger uncertainty of the influent temperature, up to +0.6°C, propagate to the model results as it was used as input variable and thereby added to the overall model uncertainty. A general uncertainty analysis of all major input variables was not conducted in this study (Belia et al. 2009).
The net average temperature increase over the plant (from influent to effluent) for the full year was 0.78°C ( Figure 5(a)). However, seasonal variations occur (Figure 3). During the warmer period of the year, the change was largerdue to solar radiationwith a net increase in temperature ( Figure 5(b)). For example, the temperature increased by 1.8, 1.4 and 1.2°C for June, July, August 2019, respectively. During winter, the net changes were small, À0.13, À0.04 and À0.05°C for November, December 2019 and January 2020, respectively. With minimal solar radiation in winter (Figure 2(c)), the energy contribution from the biological processes at this time of the year was not even compensating for the energy loss due to convection, evaporation, aeration etc. (Figure 6(b)). Studying the sub-processes of the plant individually, the ASU was the one with the greatest impact on temperature (Figure 5(a)). Excluding the ASU, the average temperature change was 0.35°C, or 45% of the total change. The temperature change in the primary treatment steps (grit chamber, pre-aeration and primary clarifiers) comprise 6% of the total. This is not insignificant and shows the importance of evaluating the temperature impacts plant-wide for accurate simulations of biological processes, for example in the ASU. Figure 6(a) shows that the heat flux components with the greatest impact on average was f sr and f bp (positive) and f ar and f c (negative). The process with the clearest seasonal variation was f sr (b), which follows the measured solar radiation (c). Also f c and f e showed variations from month to month with correlation to the outside temperature (a). The expression for f c in the model reflects energy transfer for forced convection with mixed boundary layer conditions (2:5 Á 10 5 , Re , 1:5 Á 10 7 with less than 10% of values , 5 Á 10 5 ). Whilst this provided the best fit to data of the alternatives tested (results not shown) for the presented case study, ambient conditions (e.g. wind speeds) in other cases might require altering this expression. For the energy contributions by biological reactions (f bp ) the degradation of COD is dominating while nitrification and denitrification have smaller and similar impacts. During model development, the impact of rain was evaluated. While storms heavily impact the influent flow and temperature (captured in influent data to the model) the rain over the uncovered tanks did not have a significant effect on temperature (data not shown).  Table 2 and Figure 7. The temperature increase over the plant was larger at lower influent temperatures, 1.0°C in Scenario 3 compared to 0.8°C for the default case. Considering the temperature in the ASU, which is critical for biological reaction rates, most of the temperature increase was attained there. When lowering the influent temperature 3°C the temperature at the ASU effluent was reduced 2.8°C.
In the ASU, nitrification declines as temperatures drop. The annual average secondary effluent ammonia concentration increases from 4.9 mg/l in the default scenario to 7.5 mg/l in Scenario 3. However, the impact was not linear as temperatures dropped. The increase from the default case to Scenario 1 was 0.65 mg/l and from Scenario 2 to 3, it was 1.1 mg/l. The seasonal variation was also strong (Figure 7(a)). While effluent ammonia concentrations were about the same for the three scenarios during summer, the increase was significant during winter.
The plant effluent shows a similar behaviour as the secondary effluent. The tertiary treatment with post nitrification and denitrification reduced some of the remaining ammonia and nitrate. In the plant effluent, the increase in ammonia was Figure 3 | Measured influent and effluent temperature for the treatment plant along with simulated effluent temperature for the whole validation period (a) and selected periods for summer (b) and winter (c). In Figure 3(a), the blue shaded area marks the selected periods and the red shaded area marks the period of partly corrupt data for the measured effluent temperature. In Figure 3(b) and 3(c) the red shaded area between dashed grey lines mark the measurement uncertainty for the measured effluent temperature.    Water Science & Technology Vol 00 No 0, 10 Uncorrected Proof limited for Scenario 1, small for Scenario 2 but around 1 mg/l for Scenario 3. The post-nitrification step manages a slightly increased load of NH 4 -N but not fully for Scenario 3. The post-denitrification process (with nitrate control adjusting the ethanol dosage) was able to balance the nitrate levels in the effluent. The average total nitrogen in the effluent increased up to 9.1 mg/l (Scenario 3) compared to 7.7 mg/l for the default case while the ethanol addition increased by 8%. The annual average of 9.1 mg/l TN was below the plant effluent standard which is 10 mg/l for TN. However, the dynamic simulation (Figure 7(b) and 7(c)) revealed high peaks in both NH 4 -N and TN during winter. This would not be a desired process performance at any WWTP. For Scenario 1 and mostly for Scenario 2, the TN stayed below 10 mg/l also during winter.
The different temperatures and their impacts on the microbial nitrogen conversion processes in the ASU and MBBR led to changes in plant operation ( Table 2). The plant controllers compensated for the changes in the processes. Primarily, the aeration in the ASU was reduced as less oxygen was consumed for nitrification and oxygen was more easily dissolved at the lower temperatures. Aeration was reduced by 3, 7 and 11% for Scenario 1, 2 and 3, respectively, compared to the default case. Furthermore, the carbon addition control to the post-denitrification reactor compensated for the changes in nitrate load in the three scenarios. The carbon consumption increased for all scenarios with reduced temperature as more nitrogen was pushed to the tertiary treatment from the ASU. At the most, 11% more ethanol was dosed in Scenario 2 compared to the default case. However, the trend was not distinct due to the sequential impacts on NH 4 -N and NO 3 -N concentrations from the ASU and post nitrification described above. As energy for aeration was a relatively larger cost than external carbon source (5,940 compared to 4,350 cost units per day for the default case using the operational cost index equations from BSM2), the savings on aeration exceeded the increased expenditure for denitrification. Figure 7 | Dynamic simulation results (daily averages) for ammonia in the activated sludge unit (ASU, a.) and ammonia and total nitrogen in the plant effluent (b. and c.) for the three scenarios reducing the influent temperature with 1, 2 and 3°C compared to the default case.