Resources

Start-up phaseTo avoid the toxic effect of higher nitrite concentration on microorganisms30, the influent NH4+–N and NO2−–N during the first 66 days was maintained at 50 and 60 mg/L, respectively. The NH4+–N removal efficiency in each reactor reached up to 60% quickly with a certain amount of NO3−–N in the effluent. On the 5th day, the NH4+–N removal from reactors R1–R4 was 32.1, 36.9, 40.4, and 39.1 mg/L (Fig. 1a), respectively cumulative concentration of NO3−–N was 23.7, 33.3, 39.9, and 31.1 mg/L in respective reactors (Fig. 1c). The higher concentration of NO3−–N might be due to the proliferation of nitrifying bacteria or lower denitrifying bacteria activity. It has also been reported that other bacteria, despite nitrite-oxidizing bacteria (NOB), carrying nxr gene, can contribute to nitrate production from nitrite in an anammox reactor31. Besides, concentrations of NO2−–N in each reactor’s effluent showed a different degree of an increasing trend. On the 11th day, the NO2−–N concentration in the effluent of R1–R4 reached 60.5, 52.6, 62.8, and 54.1 mg/L, respectively (Fig. 1b).Fig. 1: Nitrogen removal performance of four reactors during the start-up period.a NH4+–N concentration (mg/L), b NO2−–N concentration (mg/L), c NO3−–N concentration (mg/L).Full size imageFurthermore, the DO of each reactor was also monitored, and it was found that DO values of four reactors in the first 20 days were varied from 0.20 to 0.50 mg/L. The oxygen half-saturation constant of AOB has been reported in the range of 0.2–0.4 mg/L, and the growth rate of AOB can reach 0.5–1.0 per day32. Thus, AOB and archaea, and some NOB adapt to the lower DO concentration and could grow33,34,35. So, the higher NO3−–N values in the start-up phase can be justified by the growth of nitrifying bacteria and some other bacteria carrying the nxr gene. However, as the incubation time increased, some bacteria which may not be adapted to the conditions might be degraded. The denitrifying bacteria could propagate on the carbon source provided by degraded biomass12. The use of NOx−–N by denitrifying bacteria contributed to the nitrogen removal rate (Figs. 1 and 2c) and can be called denitrifying propagation phase. Due to unfavorable environmental conditions, the degraded bacteria could have provided a carbon source for denitrification bacteria’s growth13,14. It has been seen from Fig. 2c that in the initial, nitrogen removal efficiency showed a decreasing trend, which is caused by consumption of endogenous organic matter, which slowly consumed and excluded the dominancy of the denitrification process. The conditions in the reactors are beneficial to anammox growth and eventually led to the washout of denitrifiers. The high TN removal in the latter phase was attributed to the enrichment of anammox abundance.Fig. 2: Nitrogen removal performance of the four reactors in the stabilization period.a NH4+–N concentration (mg/L), b NO2−–N concentration (mg/L), c removal efficiency (%) of total nitrogen.Full size imageFrom the 29th day, the effluent concentrations of NH4+–N and NO2−–N in R3 declined sharply, and the NO2−–N and NH4+–N removal ratio increased from 0.52 to 1.24, close to the theoretical ratio of 1.32, which indicates a rapid increase in the anammox bacterial activity5,36. The stoichiometric ratios of R1, R2, and R4 reached close to theoretical value on days 50, 50, and 41, respectively. The activity enhancing phase started around 22 days earlier in R3 when compared with R1, which clearly showed the positive effects of MF to short the start-up period of anammox. Though R4 showed a more extended start-up period than R3, it was still shorter when compared with R1. It was observed that R2 with the only nZVI did not show a clear difference from the control reactor. On the 34th day, the NH4+–N and NO2−–N removal rates of R3 were 95.5% and 88.1%, respectively, which maintained at 85% from days 41 to 60. Compared with the control, the start-up time of R3 was shortened from 60 days to 34 days under the action of MF, and the efficiency was increased by 43.3% in the start-up period. The start-up time of the reactor under the combined effects of nZVI and MF was ~50 days, and the efficiency increased by 16.7%. On the other hand, the start-up time of anammox in R1 and R2 showed no noticeable difference.Reactors’ performance under increased nitrogen loadingAfter the successful start-up of all reactors, the influent NH4+–N and NO2−–N concentrations were increased progressively to observe the reactors’ stability under different loading rates. As shown in Fig. 2, on day 69, the influent concentrations of NH4+–N and NO2−–N were increased from 50 to 100 mg/L and 60 to 120 mg/L (nitrogen loading rate was 2.64 kg/L/d), respectively. The response of R2 against increased N loading was almost comparable to R1. The concentration of effluent NO2−–N in R1–R3 was slightly increased. The removal performance of R4 was somewhat better than R1. Similarly, the removal of NH4+–N and NO2−–N in R4 was better than R3 (Fig. 2a, b). This implies that the increase in nitrogen loading does not have a deterioration impact on the operation of the four reactors after the successful build-up of anammox activity.The influent concentration of NH4+–N and NO2−–N was again amplified to 200 and 240 mg/L (nitrogen loading rate was 5.28 kg/m3/d) on the 91st day. After the second increment in N loading, the effluent NH4+–N concentration in R1 and R4 was increased to 69.2 and 56.3 mg/L, respectively, and the effluent NO2−–N concentration was increased to 55.5 and 37.4 mg/L in R1 and R4, respectively. The increasing trend in effluent NO2−–N concentration in R1 was observed till the 106th day, and the highest value of effluent NO2−–N concentration was 134.8 mg/L which decreased the removal rate in R1 about 42.4%. Likewise, considerable fluctuations were also observed in the removal rate of NH4+–N, where more variation was noted in R1. A little higher nitrogen removal efficiency of R3 on day 91 (80%) and on day 106 (71.18%) were witnessed compared to R4 (67% and 66%, respectively). Similarly, a little high nitrogen removal efficiency was registered in R2 on days 88–91 and on days 97–100 as compared to R4. These variations might be caused due to the calculation error. The presence of MF and nZVI decreased the influence of substrate shock on the anammox performance. Although the removal rates of NH4+–N and NO2−–N in R4 were also decreased, but the fluctuation was slighter than R1 and R3. On the 115th day, the removal rates of NH4+–N and NO2−–N in R4 were recovered to over 80% and gradually reached over 90% afterward. Comparable results have also been observed by Wang et al.14 in ABBR. However, Chen et al.6 reported different results, which showed anammox was failed to adapt to the higher nitrogen loading rate. In addition, the increase in nitrogen loading also showed a negative impact on the nitrogen removal performance of R3. As shown in Fig. 2c, total nitrogen removal efficiency (about 80%) of R4 was always better than other reactors from the 90th to 180th day.Fig. 3: Nitrogen mass balance of R4 on day 166.This mass balance is drawn on the base of anammox reaction stoichiometry. The black arrows showed the amount of nitrogen to effluent, the red lines indicated anammox process, the green line indicated partial nitrification process, and the pink arrows indicate partial denitrification or partial-DNRA process.Full size imageNitrogen mass balance (Fig. 3) on day 166 proved the coupling of nitrogen removal bacteria on the basis of reported anammox stoichiometric values of NO2−–N/NH4+–N (1.32) and NO3−–N/NH4+–N (0.26). So, if there is only an anammox process, ammonium and nitrite consumption should be around 196 and 258 mg/L, respectively, to meet the reported stoichiometry of the anammox process and effluent nitrate concentration should be approximately 51 mg/L. However, the stoichiometric values of NO2−–N/NH4+–N and NO3−–N/NH4+–N obtained in this study were around 1.15 and 0.19, respectively, in R4 (coupled effects of nZVI and MF) on day 166. This means, 227 mg/L of NO2−–N and 196 mg/L of NH4+–N consumed by anammox, and the remaining amount is removed by other nitrogen cycle bacteria. The NH4+, NO2−, and NO3− attributed to various groups of microorganisms by keeping R1 (1.15) and R2 (0.19) values by using Eqs. 4 and 5. Overall, the proposed nitrogen mass balance shows a clear coupling of anammox with other nitrogen cycle bacteriaEffects of nZVI and MF on the functional genesSo as to study the combined effect of nZVI and MF on functional genes of anammox bacteria, the 16S rRNA gene copy numbers of anammox and other nitrogen cycle bacteria at different stages of the experiment were determined by the qPCR technique. The copy numbers of hzo (hydrazine oxidoreductase) functional gene were also recorded.The anammox 16S rRNA copy number was increased gradually with time (Fig. 4a). Initially, the copy number of anammox 16S rRNA was 1.01 × 106 copies/ng DNA. R1 had the lowest anammox copy number among the four reactors on the 60th day, 120th, and 180th day. The highest copy number of anammox 16S rRNA was recorded in R4 (2.13 × 106), which was followed by R3 (2.10 × 106) and R2 (1.70 × 106) on the 60th day. Nevertheless, R4 and R3 have almost similar gene copy number on the 60th day, which made us propose that MF positively influence the anammox activity, which reduced the start-up period of anammox. It has been reported that bacterial activity is suppressed in the presence of incompatible nZVI concentrations37,38,39. The adaptive capacity of microbes depends on the nZVI concentration40. The inhibition in anammox activity at higher nZVI concentration (3 g/L in this experiment) is also supported by a comparable study41. Interestingly, at day 180, the copy numbers of anammox 16S rRNA in R2 (7.17 × 106 copies/ng of DNA), R3 (6.98 × 106 copies/ng of DNA), and R4 (7.8 × 106 copies/ng of DNA) were significantly higher than R1 (5.56 × 106 copies/ng of DNA). It has also been claimed that the optimum concentration of nZVI can improve the proliferation of anammox cells after adaptation41. It is reported that nZVI lost the reactivity after 3 and 60 days in the presence and absence of oxygen in the wastewater treatment system, respectively42. It is an accepted fact the nZVI released Fe2+ and H2 in the solution43, which can be stored by anammox. Anammox can store iron ions for future haem synthesis and haem-containing enzymes involved in the electron transport chain44. Further, according to the previous study, the addition of nZVI enhanced the abundance of anammox bacteria22,23,45. The above explanation supports the conclusion about the higher anammox gene copy number in R2 in the later stage. Compared with R1, the percent increase of anammox gene copy numbers under a higher nitrogen loading rate were 29.0%, 25.5%, and 40.3% in R2–R4, respectively, on day 180.Fig. 4: The qPCR results of anammox enrichment process.a Anammox 16S rRNA copy numbers at different days (0, 60, 120, and 180) and b functional gene copy numbers of AOB, denitrifying bacterial, and anammox on day 180. Data indicate average, and error bars represent standard deviation of the results from three independent samplings, each tested in triplicate.Full size imageFurther investigation was done to analyze the functional gene hzo of anammox and other nitrogen-cycle related bacteria in the reactors (Fig. 4b) on day 180. The co-existence of AOB and denitrifying bacteria with anammox is reported to improve nitrogen removal efficiency46. The copy number of hzo gene in R1–R4 were 2.39 × 105 and 3.87 × 105, 3.29 × 105, and 4.01 × 105 copies/ng DNA, respectively. The coupled effects of nZVI and MF on anammox functional gene copy number were obvious, which might be responsible for the higher nitrogen removal under increasing nitrogen loading. The contents of amoA, nirK, and nirS in R3 were less, different from other reactors. The difference in amoA, nirK, and nirS gene copy numbers implies that microbes behave differently under the influence of MF. The optimum range of MF intensity is dissimilar for different microorganisms28,47, which can be considered a possible reason for the difference in gene contents among all four reactors. The lower nirK and nirS gene contents in R3 indicated that MF of this intensity (65 ± 10 mT) might have inhibitory effects on denitrifying bacteria. Furthermore, higher nirK gene contents in R4 (nZVI + MF) originate a premise that nZVI might have neutralized the adverse effects of MF on denitrifying bacteria. As it can be seen from Fig. 4b, the denitrifying functional genes nirK and nirS were higher in nZVI reactor. Further, the lower amoA gene content signifies anammox and denitrifying bacteria as major contributors in nitrogen removal.Effects of nZVI and MF on microbial communitiesIt is an established fact that community structure and abundance considerably affect the stability and performance of the anammox process. Therefore, the community composition of the four reactors was analyzed after 180 days by high-throughput sequencing.The sludge samples from all four reactors were collected at the end of the experiment, and operational taxonomic units (OTUs) and four indices (ACE, Chao1, Simpson, and Shannon) of each sample were calculated (Table 1). The OTUs varied from 1599 to 1889 for different treatments. The ACE and Chao1 indices reflect the community’s richness, while Simpson and Shannon’s indices reflect the community’s diversity48. Generally speaking, the greater the Chao1 or ACE index, the higher the abundance of the community. The Shannon diversity index comprehensively considers the richness and uniformity of the community. The higher the Shannon index value, the higher the diversity of the community. The Simpson index is also one of the commonly used indexes for evaluating community diversity. The higher the Simpson index value, the higher the community diversity. The value of ACE (1600) and Chao1 (1600) indices of R1 was the minimum in four reactors reflecting the lowest community richness in R1. On the other hand, the maximum community richness was observed in R3 (MF reactor) as reflected by the highest value of ACE (2180) and Chao1 (2145) indices followed by R2 (nZVI reactor) and R4 (nZVI + MF reactor). In short, the effect of MF on community richness was more pronounced than the effect of nZVI and the combination of nZVI and MF. The Simpson indices of all samples were not different. However, the Shannon indices of R4 were slightly higher than other samples, which showed a little higher diversity of community than other reactors.Table 1 The OTU numbers and bacterial diversity indices of sludge samples.Full size tableThe effects of nZVI and MF treatments on community structure distribution of sludge samples at different classification levels are presented in Fig. 5. The Chloroflexi, Proteobacteria, Cholorobi, and Planctomycetes were the abundant phyla in all samples with little variation among different treatments (Fig. 5a). The presence of phyla such as Chloroflexi, Proteobacteria, and Cholorobi in an anammox reactor was also reported previously49. Relative abundance of Chloroflexi in R1–R4 was 24.3%, 29.5%, 23.4% and 30.0%, respectively. The phylum Chloroflexi has also been detected in the anammox reactor and single-step autotrophic nitrogen removal system50,51. Some genera’s proposed role belongs to Chloroflexi in the anammox reactor is to consume the dead organic material and avoid their accumulation52. The relative abundance of Proteobacteria in R1–R4 was 25.2%, 28.5%, 22.5%, and 27.5%, respectively. According to the previous literature, mostly nitrifying and denitrifying bacteria (Thauera, Denitratisoma, and Geobacter) belong to the phylum Proteobacteria and could use NO2+–N and NH4+–N for their metabolism and proliferation activities4. As far as the relative abundance of phylum Planctomycetes in different reactors is concerned, no significant difference in R1 (8.4%), R3 (8.4%), and R4 (8.5%) was noted in this regard (chi test, p  More

The global databaseAs of 10 December 2019 a total of 9503 HAEDAT events had been entered from across the globe, comprising 48% seafood biotoxin, 43% high phytoplankton counts and/or water discolorations causing a socio-economic impact, 7% mass animal or plant mortalities and 2% others (including foam and mucilage production). In a number of HAEDAT records, a single incident was categorised into multiple event types, such as both water discoloration and high phytoplankton count (11% were multiple event types). Among all events linked to seafood toxin syndromes, Paralytic Shellfish Toxins (PST) accounted for 35%, Diarrhetic Shellfish Toxins (DST) 30%, Ciguatera Poisoning (CP) and marine and brackish water cyanobacterial toxins each 9%, Amnesic Shellfish Toxins (AST) 7% and others 10% (including Neurotoxic Shellfish Toxins (NST), Azaspiracid Shellfish Toxins (AZA) and toxic aerosols).Different OBIS regions (Fig. 1A) contained varying numbers of HAEDAT reports, with the largest number of records available for Europe, followed in order by North Asia, Mediterranean, the east and west coast of North America, the Caribbean, Pacific/Oceania, South East Asia and more limited data sets for South America, and Australia/New Zealand (Supplementary Table 1).All geographic regions were impacted by multiple HAB types, but in varying proportions (Fig. 1B). High phytoplankton density problems comprised more than 50% of regional HAEDAT records in the Caribbean, Benguela, Mediterranean Sea, North and South East Asia, while seafood toxins and fish kill impacts dominated in all other regions. The productive Benguela upwelling system is prone to mass benthic mortalities linked to high algal biomass and low oxygen19. Among toxin-related impacts, those caused by PST prevailed in North America, the Caribbean, South America, South East Asia and North Asia, whereas DST were the most frequently recorded in Europe and the Mediterranean (Fig. 1C), and are an emerging threat in the USA. NST were confined to Florida (US), with a single outbreak also reported from New Zealand. Human poisonings from Ciguatera were prominent in the tropical Pacific, the Indian Ocean, Australia and the Caribbean. While HAEDAT records of PST, DST and AST mostly relate to the incidence of phycotoxins in seafood, causing closures of shellfish harvesting areas but rarely being associated with human poisonings, records for Ciguatera (CP) refer almost exclusively to human poisonings reported by medical practitioners. HAB events related to marine cyanotoxins were mainly reported from Australia/New Zealand (ANZ), the Indian Ocean (IND) and the Baltic Sea area in NW Europe (EUR).Temporal trends in HAB eventsThe number of HAEDAT events reported for each region per year generally showed increases (Fig. 2A). Specifically, eight of the nine regions showed increases of which six were statistically significant (ECA, CCA, WCA, SEA, MED, EUR; Supplementary Table 1). The meta-analysis of the relationship between HAEDAT events and Year, plotted for each region with the significance level and the confidence intervals adjusted by the effective degrees of freedom (reduced because of autocorrelation) is shown in Fig. 2B. Most of the correlations for the individual regions overlapped zero (i.e. they were not significant), but the overall global total meta-analytic correlation was significant (r = 0.37, z = 2.97, p = 0.003) because it combined the information from each individual region and thus provided more statistical power. This suggests that the number of HAEDAT events is increasing over time. The meta-analysis of the relationship between number of geographic grids with one or more HAB events and Year had similar results but weaker statistical significance (r = 0.27, z = 2.16, p = 0.031, Fig. 2C).While the number of geographic grids with HAB events is less prone to inconsistencies in what constitutes a HAB event and less affected by sampling effort, to try and adjust more specifically for sampling effort we used OBIS data on microalgae sampling. OBIS data generally showed an increase in sampling effort, although the SAM and NAS regions did not follow this trend (Fig. 3A). Once HAEDAT events were adjusted relative to OBIS observations, there were contrasting trends in standardised HAB events over time, with four regions (SAM, WCA, ANZ and NAS) changing direction compared to the unadjusted HAEDAT data (Fig. 3B vs. Fig. 3A). The meta-analysis of the standardised HAB events showed five of the nine regions with a substantially increasing trend (two flat, and two declining), but there was no statistically significant trend overall when all regions were combined (r = 0.35, z = 1.33, p = 0.18), (Fig. 3C). This implies there is insufficient evidence to conclude that HABs are increasing across all the regions analysed, but it is clear there are contrasting trends in individual regions.In the period studied, aquaculture production increased 16-fold from a global total 11.35M tonnes in 1985 up to 178.5M tonnes in 2018, with the largest increases occurring SEA and SAM + CCA and with North America (ECA + WCA) and EUR stabilising (Fig. 4A). The number of HAEDAT events over time was significantly correlated with aquaculture production, with all regions with suitable data exhibiting more HAEDAT events as aquaculture expanded, with a strongly significant relationship overall (r = 0.43, z = 3.59, p = 0.0003; Fig. 4B bottom; Supplementary Table 2).Fig. 4: Changes in different geographic regions of aquaculture production in the period 1985 to 2018, and meta-analysis of HAEDAT events over time against aquaculture.A Changes in 1985 to 2018 in five regions (ECA + WCA; SAM + CCA; ANZ; SEA; EUR) of tonnage of Aquaculture Production of fish, molluscs, crustaceans and aquatic plants; and B. Meta-analysis of HAEDAT events over time vs. Aquaculture. The overall number of HAEDAT events over time was significantly correlated with aquaculture production (bottom). Weighted mean correlations (filled circles) are shown with 99% confidence limits (bars) in (B).Full size imageSelected HAB case studiesFurther exploring the influence of monitoring efforts, Fig. 5 depicts a 4× fold increase of positive global records between 1985 and 2018 of the main causative organisms of Diarrhetic Shellfish Poisoning (DSP; 84,392 OBIS records of the dinoflagellate genus Dinophysis; Fig. 5A), a 7× fold increase of global observations of the main causative organisms of Amnesic Shellfish Poisoning (ASP; 128,282 records of the diatom genus Pseudo-nitzschia; Fig. 5B) and 6× fold increase of global observations of one of the causative organisms of Paralytic Shellfish Poisoning (PSP; 9887 records of the dinoflagellate genus Alexandrium; Fig. 5C). It should be noted that records for Dinophysis, Pseudo-nitzschia and Alexandrium may also include non-toxic species or strains. In all three cases the clear increase in the number of observations through time of causative organisms is paralleled by the increase of HAEDAT records of the associated toxin syndromes (Fig. 5D–F) which in the case of PSP are contributed also by other species, namely the tropical Pyrodinium bahamense and widespread Gymnodinium catenatum. The occurrence of the causative toxigenic HAB species is not always an accurate predictor for the incidence of human shellfish poisonings (indicated by the shellfish icons in Fig. 5G–I). This reflects the efficiency of the food safety risk management strategies implemented in many of the affected countries. Globally, some 11,000 cases of nonfatal events related to DSP were reported, mostly from Europe, South America and Japan (Fig. 5G). These events mainly include closures of shellfish harvesting areas due to observed levels of DST above regulatory limits to protect human health. It is noted that, despite the widespread distribution of Pseudo-nitzschia species (Fig. 5H), there have been no human fatalities from Amnesic Shellfish Poisoning since the original 1987 incident in Prince Edward Island, Canada (150 illnesses with 3 fatalities), even though associated mortalities of marine mammals of high conservation value are of increasing concern in western North America, including in the climate hotspot of Arctic Alaska36. AST has also been associated with marine mammal calf mortalities in Argentina37. Of the global total of 3800 human Paralytic Shellfish Poisonings during the 1985–2018 period (Fig. 5I), the largest number occurred in the Philippines, a country strongly dependent on aquaculture for human food protein, with 2555 poisonings recorded between 1983 and 2013 of which 165 were fatalities23,38, predominantly caused by highly toxic Pyrodinium bahamense. Because of the increased use of molecular detection methods our knowledge on the global distribution of ciguatera- causing organisms, selected species of the benthic dinoflagellates Gambierdiscus and Fukuyoa, has increased considerably (Fig. 6A). The database on the presence of ciguatoxins in fish (Fig. 6B) is still limited because of the complexity of the chemical analysis used to confirm the presence of ciguatoxins. Exploring trends of human CP, in Hawaii poisonings have been decreasing, in French Polynesia and the Caribbean numbers remained stable, whereas CP is a new phenomenon in the Canary Islands (Fig. 6C). Globally CP affects 10,000–50,000 people annually but fatalities are rare20.Fig. 5: Increases between 1985 and 2018 of global observations of the causative organisms, HAEDAT toxic events, and distributions of the toxin syndromes Diarrhetic, Amnesic and Paralytic Shellfish Poisoning.A–C Total number of global observations from OBIS of causative microalgal organisms of Dinophysis spp., Pseudo-nitzschia spp. and Alexandrium spp.; D–F The number of records of HAEDAT Toxic Events of DST, AST and PST.; G–I Global distribution maps (as red dots, from OBIS) as of 2018 of Dinophysis spp. (DSP), Pseudo-nitzschia spp. (ASP) and Alexandrium, Pyrodinium, Gymnodinium catenatum (PSP). The locations of toxic events resulting in human poisonings are indicated by the size of the shellfish icons. The first number shows number of poisonings, the second number indicates fatalities. For ASP, 150/3 signals 150 clinical cases with three fatalities. No human fatalities have ever occurred from DSP.Full size imageFig. 6: Known global distribution in 2018 of the causative dinoflagellate genera, ciguatoxins in fish, and trends of human ciguatera poisonings in selected geographic regions.A Distribution of the dinoflagellate genera Gambierdiscus and Fukuyoa (blue and orange dots); B Ciguatoxins in fish (red) and shellfish (orange); and C Trends between 2000 and 2018 in human ciguatera poisonings in Hawaii, French Polynesia, Canary Islands, the Caribbean (light green) and Mexico (dark green). Adapted from Chinain et al.20.Full size imageAquacultured finfish mortalities caused by the taxonomically unrelated microalgal genera Chattonella, Pseudochattonella, Heterosigma, Karenia, Karlodinium, Margalefidinium (Cochlodinium) and Prymnesium/ Chrysochromulina globally account for much greater economic damage than HABs contaminating seafood39. While most shellfish toxins have now been well characterised and are effectively monitored and regulated, finfish held captive in intensive aquaculture operations continue to be vulnerable to HABs (USD71M loss in Japan in 1972, USD70M in Korea in 1995, USD290M in China in 2012, USD100M in Norway in 201919,20,21,22,23,24,25,26,27,28,29,30,31,32), even though the causative ichthyotoxins usually are of no human health significance. The 2016 Chilean salmon mortality that caused a record USD800M loss led to major social unrest40. Again, the incidence of fish-killing HAB species is not an accurate predictor of economic losses. For example, Heterosigma blooms occur both on the west and east coast of North America, but fish mortalities are mostly confined to the west coast29. In large part, this reflects locations where blooms occur relative to the location and size of the aquaculture operations. The dinoflagellate Karlodinium australe never caused any problems in its Australian lagoon type locality41 but in 2014 killed 50,000 caged fish in Malaysia and is now also known from Japan and the Philippines42. In the wild, finfish can swim away from bloom areas, hence aquaculture finfish mortality is largely a human-generated problem. More

Monday is sample-collection day in Cape Town, South Africa, and Aqeelah Benjamin is halfway through her shift. At the Green Point wastewater-treatment plant, under the Atlantic coast promenade, she fills a 500-millilitre bottle from a tap of untreated water. She wipes the bottle’s exterior, cleans it with a spritz of ethanol and stores it on ice.It’s one of nine samples that Benjamin will collect today, each from a different facility. Later, she’ll drop them off at the South African Medical Research Council (SAMRC) laboratory, where they will be tested for the presence of SARS‑CoV-2, the virus that causes COVID‑19. But that’s just a fraction of the samples’ potential — waste water contains a wealth of information about public health, and scientists are only just starting to tap into its potential.Wastewater testing can provide an unbiased snapshot of community health: whatever access they have to the health-care system, everyone has to go to the toilet. And as the effluent makes its way to treatment facilities, researchers can test it to determine what pathogens might be present. For more than 40 years, researchers have used waste water to monitor the spread of poliovirus. Norovirus, influenza, hepatitis and measles viruses can also be found in waste water. Now, a growing number of countries are using waste water to monitor the spread of COVID-19. The memorably named COVIDPoops19 Dashboard, run by researchers at the University of California, Merced, lists more than 2,200 monitoring sites in 54 countries (see go.nature.com/3fjfcjt).South Africa is one of a handful of countries rolling out the technology nationwide. But the process isn’t easy: researchers are struggling to overcome logistical hurdles and extend the techniques to the large part of the population that has no sewerage infrastructure.
The myriad ways sewage surveillance is helping fight COVID around the world
The SAMRC runs a research programme across four of South Africa’s nine provinces — the Western Cape, home of Cape Town; the rural Eastern Cape and Limpopo; and Gauteng, which includes South Africa’s largest city, Johannesburg, and its administrative capital, Pretoria. After showing that it could detect SARS-CoV-2 in waste water in 5 treatment plants last June, the SAMRC extended the testing to another 19 plants to work out the logistics of scaling up the work.The exercise highlighted problems specific to operating in South Africa. Difficulties collecting samples from remote sites can slow down the diagnostic process, for instance. And rolling electricity blackouts, known locally as ‘load shedding’, can hinder the operation of the equipment that samples waste water throughout the day. These machines, known as composite samplers, are also prone to theft. Considering this, and the cost of the samplers, South African researchers tend to use ‘grab samples’, such as Benjamin’s. These samples are generally considered less representative than are those from composite samplers, because they represent just a single snapshot in time. But a study by the South African Water Research Commission (WRC) has found little difference in the effectiveness of the two approaches in detecting SARS-CoV-2 (see go.nature.com/3v1mpm4).On Monday afternoon, after Benjamin and two other collectors deliver their samples to the SAMRC lab in the northern Cape Town suburb of Parow, I meet Rabia Johnson, deputy director at the SAMRC’s Biomedical Research and Innovation Platform. The lab specializes in molecular biology and cell-culture systems, and has been testing waste water for SARS-CoV-2 for nearly a year. “I think we’ve got the most comprehensive longitudinal database [in South Africa], from the first wave through the second wave,” Johnson says.

Researchers process samples for testing at their lab in Tshwane.Credit: Delwyn Verasamy/Mail & Guardian

In the lab upstairs from Johnson’s office, the researchers concentrate the samples in a centrifuge and then extract any viral RNA using a kit from the molecular-reagents company Qiagen in Hilden, Germany. The kit is optimized for extracting RNA from soil rather than water, but researchers at Yale University in New Haven, Connecticut, have shown that it is better at handling the unwanted organic materials found in wastewater than are conventional techniques (J. Peccia et al. Preprint at medRxiv https://doi.org/gc9k; 2020). The team then moves the extracted RNA to a ‘clean’ room to test for SARS-CoV-2 to avoid the risk of contamination. The researchers use a technique called real-time quantitative polymerase chain reaction (RT-qPCR) to quantify the amount of RNA that encodes the viral nucleocapsid protein. Other viral sequences are added in to assess performance. And positive controls are added for two key variants: 501Y.V2, first identified in South Africa, and B.1.1.7, detected in the United Kingdom. A fluorescent signal indicates that the nucleocapsid RNA is present.Finally, Johnson cleans up the data and sends them to the SAMRC’s Environment and Health Resource Unit. Researchers there upload the data to the SAMRC dashboard, a public resource launched in November 2020 that plots virus spread on a map (see go.nature.com/3ukn74u). Around 700 people per week access the service, according to Renée Street, a senior scientist at the unit.Early warning systemBecause wastewater testing can capture the presence of the virus 7–14 days before the onset of symptoms, it can provide valuable early warning of localized outbreaks. Health officials can then make sure the necessary resources, equipment and personal protective equipment are available, says Johnson.
How sewage could reveal true scale of coronavirus outbreak
But that’s still theoretical. Wastewater data have not been used directly to inform decisions about control measures such as targeted lockdowns or resource allocation in South Africa, but they have been used alongside other sources of information, including case numbers and hospital admissions. “It’s still very new science,” says Street. “We’re still working out what the different signals are, and at what signal level we would need to take action.”The ability of the technology to identify hotspots is governed by the service area of the treatment plant: the wider the spread, the harder it is to pinpoint small outbreaks, and at least one facility serves more than 850,000 people. But Cape Town epidemiologist Natacha Berkowitz, says that the goal is to “localize infection down to a small unit area, like a suburb”.After the pandemic, the city plans to use waste water to regularly monitor for polio and other viruses. Although South Africa has been declared polio-free, missed vaccinations during the pandemic could result in an outbreak. “If we get a positive signal for polio, we’ll look at that specific drainage area, and perhaps do increased vaccinations or community education,” Berkowitz says.Remote areasThe advantage of wastewater testing is that it samples the population without requiring any action from individuals. It’s also cost-effective, because a single sample can be tested for multiple pathogens. And most labs can do the testing. “If you have a medical set-up for pathogen testing,” says Janet Mans, a virologist at the University of Pretoria, “you should be able to do this.”But the technology also has downsides. It’s not easy to tie a signal to a specific location or group, for instance. Furthermore, it monitors only households that are connected to the sewerage system. And some 40% of the nation’s households do not have a flush toilet attached to the sewers, according to the 2011 census. Extending testing to these areas would capture much more of the population.

A centrifuge is used to concentrate the samples.Credit: Delwyn Verasamy/Mail & Guardian

In Pretoria, a private facility known as Waterlab is working with the WRC to expand testing to unsewered communities. The idea is to build a framework to start using samples from rivers and surface waters, says Gina Pocock, Waterlab’s specialist consultant.To monitor trends over time, samples are taken from sites that are consistently contaminated with waste water. That includes rivers downstream of unsewered, informal (or unplanned) settlements and surface run-off both of ‘grey’ water from bathing and sinks around standpipes and of ‘black’ water, which pools near latrines and contains faecal matter.Such samples can be difficult to work with. Mans, who is working with Waterlab on ways to extract and test nucleic acids, says that environmental samples often contain compounds that can inhibit the amplification of nucleic acid by PCR, so addition of internal control sequences are a must. If the internal control is still inhibited after the sample has been diluted by one part to ten, that sample is considered invalid, says Mans. A target can be considered negative only if the internal control amplifies at that dilution.Equally difficult is the analysis, especially determining how many people the sample might represent. The researchers have to use overall trends in other parameters as proxies for the number of people. At Waterlab, Pocock says, researchers use bacterial density “to get an indication of the faecal load in the water”.We don’t know “how many people flushed their toilet this morning”, Pocock says. “It’s not a definite science, where you can say this is our viral load and X amount of people in this community are sick. And with the rivers, even less so. So, we look at trends.”Pros and consThose trends should help in assessing the effectiveness of South Africa’s response to COVID-19.Rolf Halden, director of the Biodesign Center for Environmental Health Engineering at Arizona State University in Tempe, has been testing waste water for nearly 20 years. Last year, he took part in a study to look at the feasibility of mass surveillance, testing the waste water of 36 million people in 100 US cities twice a week for 8 weeks for SARS-CoV-2. His goal is to scale that up to one billion people globally. He and his team found that it was possible to collect a lot of actionable information for very little investment, while still protecting people’s privacy.
NatureTech hub
Although obviously enthusiastic about the technology’s potential, Halden acknowledges its shortcomings. For instance, the temperature at the monitoring site matters, as does the distance that people live from it. A signal at a monitoring site could come from a single person close to the site, or from 10,000 people farther away, he explains.There are also ethical and moral considerations as the technology becomes more powerful. Aggregated data on populations are generally considered safe from a privacy perspective, because individuals cannot be identified. However, as analytical techniques advance, it might become possible to identify human DNA, prompting concerns about who should have access to both the technology and the data (D. Jacobs et al. IEEE Trans. Technol. Soc. https://doi.org/gc9m; 2021). “The moral and ethical framework has to grow, ideally before the technology is applied,” Halden says.Still, the potential benefits remain powerful motivators. Wastewater testing, Halden says, allows researchers to keep “a finger on the pulse of humanity”. More

Data for estimating drawdown areasThe calculation of drawdown areas was based on monthly time series of surface-area values for 6,818 reservoirs provided by GRSAD18. It comprises all reservoirs from the Global Reservoir and Dam dataset19 except of 45 reservoirs without reported geometric information. In accordance with ref. 1, we further removed 24 reservoirs classified as natural lakes that have been modified with water regulation structures (this includes lakes Victoria, Baikal and Ontario). The GRSAD dataset comprised entries from March 1984 to October 2015. To have a constant number of data points per year, we restricted our analysis to the period from January 1985 to December 2014.GRSAD was created by correcting the Global Surface Water dataset31 for images contaminated with clouds, cloud shadows and terrain shadows. With this correction, the number of effective images that can be used in each time series has been increased by 81% on average. Substantial improvements have been achieved for reservoirs located in regions with frequent cloud cover and high-latitude reservoirs in the Northern Hemisphere, where low illumination has previously resulted in missing area values during winter months.Calculation of drawdown areasWe calculated monthly drawdown areas for all reservoirs contained in GRSAD according to:$${rm{DA}}=left({{rm{Area}}}_{{rm{max }}}-{rm{Area}}right)/{{rm{Area}}}_{{rm{max }}}$$
(1)
where DA is the relative extent of the drawdown area for a given reservoir considering the current monthly surface area (Area) and the maximum area recorded during the period 1985–2015 (Areamax). We assumed that the maximum area of each reservoir recorded during the 30-year period is a valid representation of its nominal surface area (the area of the reservoir at maximum filling level).Complete filling of reservoirs was defined by a drawdown area smaller than 5% of Areamax. Because there is no uniform definition of ‘extreme drawdown’, we used the Cape Town water crisis 2018 as a reference21. The number of reservoirs experiencing extreme drawdown was estimated by averaging the number of reservoirs with drawdown areas exceeding 40%, 50%, 60% or 70% of Areamax at least once. To prevent initial filling of reservoirs being identified as extreme drawdown, 791 reservoirs built during the analysed period (year built ≥ 1985) were excluded from this analysis. The upper bound (70%) corresponds to the drawdown-area extent during the Cape Town water crisis 201821 (Fig. 1a). The lower bound (40%) corresponds to a reservoir capacity (storage water volume) of approximately 35%, as remained available during that water crisis, assuming an idealized, triangular reservoir shape (Extended Data Fig. 8). This was estimated according to:$$0.36=frac{{left(0.6times sqrt{2}right)}^{2}}{2}$$
(2)
For the calculation of total global drawdown area, used for the upscaling of GHG emissions, we combined data for reservoirs larger than 10 km2 with values derived from a Pareto model for smaller reservoirs. First, we estimated total reservoir surface area for nine size classes following a Pareto distribution. Subsequently, we estimated total drawdown area for each size class by multiplying the size-class-specific relative drawdown-area extent by the total reservoir surface area of each size class (Supplementary Table 3). Because the relative drawdown-area extent for reservoirs smaller than 0.001 km² is unknown and furthermore considered as being imprecise for reservoirs smaller than 10 km², we derived estimates for these size classes on the basis of four different statistical models (linear, square root, logarithmic, polynomial; Extended Data Fig. 9). Reservoirs larger than 10 km² were used to fit linear, square root and logarithmic models, whereas all available data were used for fitting a second-degree polynomial model to achieve a best representation of the available data. The four models all have a constant (linear model) or decreasing (square root, logarithmic, polynomial) slope. We have refrained from using models with increasing slopes (for example, exponential) to not overestimate the drawdown extent of small reservoirs and, thus, consider these estimates as conservative.Data analysisStatistical models to predict drawdown-area extent for each reservoir were developed using stepwise MLR. Climatic data (mean annual temperature, precipitation seasonality) for all reservoir locations were extracted from the Climatologies at High Resolution for the Earth’s Land Surface Areas climate dataset, which gives high-resolution (0.5 arcmin) climate information for global land areas over the period 1979–201332. Climate zones in the Köppen–Geiger system were determined from the high-resolution (5 arcmin) global climate map derived from long-term monthly precipitation and temperature time series representative for the period 1986–201033,34. Data on baseline water stress were extracted from Aqueduct 3.025. Baseline water stress measures the ratio of total water withdrawals to available renewable surface and groundwater supplies and is derived from high-resolution (5 arcmin) hydrological model outputs using the PCR-GLOBWB 2 model35,36.Dates were categorized into four seasons on the basis of their meteorological definition depending on hemisphere. Therefore, for the Northern Hemisphere, spring begins on 1 March, summer on 1 June, autumn on 1 September and winter on 1 December. For the Southern Hemisphere, spring begins on 1 September, summer on 1 December, autumn on 1 March and winter on 1 June.For the analyses of reservoir use types, we used the information provided in the column ‘MAIN_USE’ of the Global Reservoir and Dam dataset. Reservoirs where the main use was not specified (n = 1,554) were combined with those having MAIN_USE = ‘Other’ (n = 205).To identify the magnitude of trends in time series, we used the non-parametric Theil–Sen estimator and the Mann–Kendall test because they do not require prior assumptions of statistical distribution for the data and are resistant to outliers. The Theil–Sen estimator was used to compute the linear rate of change, and the Mann–Kendall test was used to determine the level of significance. We analysed differences between groups using the Kruskal–Wallis test and Dunn’s post hoc test. The threshold to assess statistical significance was 0.05 for all analyses, The statistical analyses were performed using R 3.4.437.Upscaling of GHG emissions and OC burialBecause the global reservoir area derived in this study differed from the area used in previous studies, we recalculated the published global estimates for both OC burial6 in and GHG emissions1 from reservoirs to allow for comparison (Extended Data Fig. 10). We fitted empirical distributions to CO2 emission data from drawdown areas (Supplementary Table 2 and Extended Data Fig. 7) as well as the published OC burial rates6 and published GHG emission data1 from water surfaces of reservoirs. For CO2 emissions from drawdown areas, we used a gamma distribution to account for non-normality of the data (Extended Data Fig. 7). For CO2 and N2O emissions from the water surface, we fitted a skewed normal distribution because of the occurrence of negative values (Extended Data Fig. 7). For CH4 emissions from the water surface, we fitted a log-normal distribution (Extended Data Fig. 7). Because the global estimate of OC burial was derived using geostatistical modelling, we fitted a gamma distribution to the published moments of OC burial rate6 (mean ± s.d. = 144 ± 75.83 gC m−2 yr−1) where the s.d. is calculated as the s.d. of the four scenarios used in that study. The final global empirical distributions for all fluxes were estimated by multiplying average emission and burial rates derived from resampling the preceding distributions times the total water surface area and drawdown area of reservoirs, resulting also from resampling their distributions after uncertainty propagation (see Treatment of uncertainty).Treatment of uncertaintyAs in all upscaling exercises, the global analysis conducted in this study is subject to substantial uncertainty. In our case, the uncertainty results from both the quantification of water surface and drawdown area of reservoirs and the estimation of global rates for GHG emission and OC burial. To comprehensively take all sources of uncertainty into account, we propagated all uncertainty throughout the whole analysis using a combination of Taylor series expansion and Monte Carlo simulations (Extended Data Fig. 10). In brief, we applied customary equations for uncertainty propagation derived from the Taylor series expansion method when propagating uncertainty of moments (for example, mean) or simple arithmetic calculations (for example, multiplication). For more-complex situations or when non-normality was conspicuous, we used Monte Carlo propagation. To obtain global estimates and standard error of water surface and drawdown area of reservoirs, both the systematic (bias) and random uncertainties of the remote-sensing-derived dataset18 as well as the uncertainty induced by our Pareto modelling for reservoirs More

1.Rowley, K. H., Cucknell, A.-C., Smith, B. D., Clark, P. F. & Morritt, D. London’s river of plastic: high levels of microplastics in the Thames water column. Sci. Total Environ. 740, 140018 (2020).CAS 
Article 

Google Scholar 
2.Hurley, R. R., Woodward, J. C. & Rothwell, J. J. Microplastic contamination of river beds significantly reduced by catchment-wide flooding. Nat. Geosci. 11, 251–257 (2018).CAS 
Article 

Google Scholar 
3.Rochman, C. M. Microplastics research—from sink to source. Science 360, 28–29 (2018).CAS 
Article 

Google Scholar 
4.Siegfried, M., Koelmans, A. A., Besseling, E. & Kroeze, C. Export of microplastics from land to sea. A modelling approach. Water Res. 127, 249–257 (2017).CAS 
Article 

Google Scholar 
5.Van Sebille, E. et al. A global inventory of small floating plastic debris. Environ. Res. Lett. 10, 124006 (2015).Article 

Google Scholar 
6.Ross, P. S. et al. Pervasive distribution of polyester fibres in the Arctic Ocean is driven by Atlantic inputs. Nat. Commun. 12, 106 (2021).CAS 
Article 

Google Scholar 
7.Napper, I. E. et al. The abundance and characteristics of microplastics in surface water in the transboundary Ganges River. Environ. Pollut. 274, 116348 (2021).CAS 
Article 

Google Scholar 
8.Long, Z. et al. Microplastic abundance, characteristics, and removal in wastewater treatment plants in a coastal city of China. Water Res. 155, 255–265 (2019).CAS 
Article 

Google Scholar 
9.López-Rojo, N. et al. Microplastics have lethal and sublethal effects on stream invertebrates and affect stream ecosystem functioning. Environ. Pollut. 259, 113898 (2020).Article 
CAS 

Google Scholar 
10.Hurley, R. R., Woodward, J. C. & Rothwell, J. J. Ingestion of microplastics by freshwater tubifex worms. Environ. Sci. Technol. 51, 12844–12851 (2017).CAS 
Article 

Google Scholar 
11.Su, L. et al. Using the Asian clam as an indicator of microplastic pollution in freshwater ecosystems. Environ. Pollut. 234, 347–355 (2018).CAS 
Article 

Google Scholar 
12.Horton, A. A., Jürgens, M. D., Lahive, E., van Bodegom, P. M. & Vijver, M. G. The influence of exposure and physiology on microplastic ingestion by the freshwater fish Rutilus rutilus (roach) in the River Thames, UK. Environ. Pollut. 236, 188–194 (2018).CAS 
Article 

Google Scholar 
13.Smiroldo, G., Balestrieri, A., Pini, E. & Tremolada, P. Anthropogenically altered trophic webs: alien catfish and microplastics in the diet of Eurasian otters. Mammal. Res. 64, 165–174 (2019).Article 

Google Scholar 
14.D’Souza, J. M., Windsor, F. M., Santillo, D. & Ormerod, S. J. Food web transfer of plastics to an apex riverine predator. Glob. Change Biol. 26, 3846–3857 (2020).Article 

Google Scholar 
15.Mateos-Cárdenas, A. et al. Rapid fragmentation of microplastics by the freshwater amphipod Gammarus duebeni (Lillj.). Sci. Rep. 10, 12799 (2020).Article 
CAS 

Google Scholar 
16.Windsor, F. M. et al. A catchment‐scale perspective of plastic pollution. Glob. Change Biol. 25, 1207–1221 (2019).Article 

Google Scholar 
17.Tang, Y. et al. Immunotoxicity and neurotoxicity of bisphenol A and microplastics alone or in combination to a bivalve species, Tegillarca granosa. Environ. Pollut. 265, 115115 (2020).CAS 
Article 

Google Scholar 
18.Woodward, J. C., Rothwell, J. J., Hurley, R. R., Li, J. & Ridley, M. Microplastics in rivers. Environ. Scientist. 29, 36–43 (2020).
Google Scholar 
19.Townsend, C. R. & Hildrew, A. G. Field experiments on the drifting, colonization and continuous redistribution of stream benthos. J. Anim. Ecol. 45, 759–772 (1976).Article 

Google Scholar 
20.Comber, S. D. W., Gardner, M. J. & Ellor, B. Seasonal variation of contaminant concentrations in wastewater treatment works effluents and river waters. Environ. Technol. 41, 2716–2730 (2019).Article 
CAS 

Google Scholar 
21.River Tame at Portwood (69027) (UK National River Flow Archive, 2020); https://nrfa.ceh.ac.uk/data/station/info/69027.html22.Upper Mersey Catchment Management Plan (National Rivers Authority, 1996).23.Sutton, R. et al. Microplastic contamination in the San Francisco Bay, California, USA. Mar. Poll. Bull. 109, 230–235 (2016).CAS 
Article 

Google Scholar 
24.Gündoğdu, S., Çevik, C., Güzel, E. & Kilercioğlu, S. Microplastics in municipal wastewater treatment plants in Turkey: a comparison of the influent and secondary effluent concentrations. Environ. Monit. Assess. 190, 626 (2018).Article 

Google Scholar 
25.Vermaire, J. C., Pomeroy, C., Herczegh, S. M., Haggart, O. & Murphy, M. Microplastic abundance and distribution in the open water and sediment of the Ottawa River, Canada, and its tributaries. Facets (Ott.) 2, 301–314 (2017).Article 

Google Scholar 
26.Lv, X. et al. Microplastics in a municipal wastewater treatment plant: fate, dynamic distribution, removal efficiencies, and control strategies. J. Clean. Prod. 225, 579–586 (2019).CAS 
Article 

Google Scholar 
27.Murphy, F., Ewins, C., Carbonnier, F. & Quinn, B. Wastewater treatment works (WwTW) as a source of microplastics in the aquatic environment. Environ. Sci. Technol. 50, 5800–5808 (2016).CAS 
Article 

Google Scholar 
28.Tibbetts, J., Krause, S., Lynch, I. & Sambrook Smith, G. H. Abundance, distribution, and drivers of microplastic contamination in urban river environments. Water 10, 1597 (2018).CAS 
Article 

Google Scholar 
29.Horton, A. A., Svendsen, C., Williams, R. J., Spurgeon, D. J. & Lahive, E. Large microplastic particles in sediments of tributaries of the River Thames, UK—abundance, sources and methods for effective quantification. Mar. Pollut. Bull. 114, 218–226 (2017).CAS 
Article 

Google Scholar 
30.Ockelford, A., Cundy, A. & Ebdon, J. E. Storm response of fluvial sedimentary microplastics. Sci. Rep. 10, 1865 (2020).CAS 
Article 

Google Scholar 
31.Blair, R. M., Waldron, S. & Gauchotte-Lindsay, C. Average daily flow of microplastics through a tertiary wastewater treatment plant over a ten-month period. Water Res. 163, 114909 (2019).CAS 
Article 

Google Scholar 
32.Curtis, E. J. C. & Harrington, D. W. The occurrence of sewage fungus in rivers in the United Kingdom. Water Res. 5, 281–290 (1971).CAS 
Article 

Google Scholar 
33.Fowles, A. E., Edgar, G. J., Hill, N., Stuart-Smith, R. D. & Kirkpatrick, J. B. An experimental assessment of impacts of pollution sources on sessile biota in a temperate urbanised estuary. Mar. Pollut. Bull. 133, 209–217 (2018).CAS 
Article 

Google Scholar 
34.Harvey, G. L. et al. Invasive crayfish as drivers of fine sediment dynamics in rivers: field and laboratory evidence. Earth Surf. Process. Landf. 39, 259–271 (2014).Article 

Google Scholar 
35.Wessel, C. C., Lockridge, G. R., Battiste, D. & Cebrian, J. Abundance and characteristics of microplastics in beach sediments: Insights into microplastic accumulation in northern Gulf of Mexico estuaries. Mar. Pollut. Bull. 109, 178–183 (2016).CAS 
Article 

Google Scholar 
36.Kowalski, N., Reichardt, A. M. & Waniek, J. J. Sinking rates of microplastics and potential implications of their alteration by physical, biological, and chemical factors. Mar. Pollut. Bull. 109, 310–319 (2016).CAS 
Article 

Google Scholar 
37.Corcoran, P. L. Benthic plastic debris in marine and freshwater environments. Environ. Sci. Process Impacts 17, 1363–1369 (2015).CAS 
Article 

Google Scholar 
38.Woodall, L. C. et al. The deep sea is a major sink for microplastic debris. R. Soc. Open Sci. 1, 140317 (2014).Article 

Google Scholar 
39.Kane, I. A. et al. Seafloor microplastic hotspots controlled by deep-sea circulation. Science 368, 1140–1145 (2020).CAS 
Article 

Google Scholar 
40.Laville, S. Water firms discharged raw sewage into English waters 400,000 times last year. The Guardian (31 March 2021); https://www.theguardian.com/environment/2021/mar/31/water-firms-discharged-raw-sewage-into-english-waters-400000-times-last-year41.Hammond, P. et al. Detection of untreated sewage discharges to watercourses using machine learning. NPJ Clean Water 4, 18 (2021).Article 

Google Scholar 
42.Acreman, M. C. & Ferguson, A. J. D. Environmental flows and the European Water Framework Directive. Freshw. Biol. 55, 32–48 (2010).Article 

Google Scholar 
43.Qiao, R., Lu, K., Deng, Y., Ren, H. & Zhang, Y. Combined effects of polystyrene microplastics and natural organic matter on the accumulation and toxicity of copper in zebrafish. Sci. Total Environ. 682, 128–137 (2019).CAS 
Article 

Google Scholar 
44.Sustainable Goal 6: Ensure Availability and Sustainable Management of Water and Sanitation for All (United Nations, 2015); https://sdgs.un.org/goals/goal6 .45.Sadoff, C. W., Borgomeo, E. & Uhlenbrook, S. Rethinking water for SDG 6. Nat. Sustain. 3, 346–347 (2020).Article 

Google Scholar 
46.Garner, G., Hannah, D. M. & Watts, G. Climate change and water in the UK: recent scientific evidence for past and future change. Prog. Phys. Geogr. Earth Environ. 41, 154–170 (2017).Article 

Google Scholar 
47.Prudhomme, C., Young, A. & Watts, G. The drying up of Britain? A national estimate of changes in seasonal river flows from 11 regional climate model simulations. Hydrol. Process. 26, 1115–1118 (2012).Article 

Google Scholar 
48.Miller, J. D. & Hutchins, M. The impacts of urbanisation and climate change on urban flooding and urban water quality: a review of the evidence concerning the United Kingdom. J. Hydrol. Reg. Stud. 12, 345–362 (2017).Article 

Google Scholar 
49.Lambert, C. P. & Walling, D. E. Measurement of channel storage of suspended sediment in a gravel-bed river. Catena 15, 65–80 (1988).Article 

Google Scholar 
50.Owens, P. N., Walling, D. E. & Leeks, G. J. L. Deposition and storage of fine-grained sediment within the main channel system of the River Tweed, Scotland. Earth Surf. Process. Landf. 24, 1061–1076 (1999).CAS 
Article 

Google Scholar 
51.Hurley, R., Rothwell, J. & Woodward, J. C. Metal contamination of bed sediments in the Irwell and upper Mersey catchments, northwest England: exploring the legacy of industry and urban growth. J. Soils Sediment. 17, 2648–2665 (2017).CAS 
Article 

Google Scholar 
52.Koelmans, A. A. et al. Microplastics in freshwaters and drinking water: critical review and assessment of data quality. Water Res. 155, 410–422 (2019).CAS 
Article 

Google Scholar 
53.Hurley, R. R., Lusher, A. L., Olsen, M. & Nizzetto, L. Validation of a method for extracting microplastics from complex, organic-rich, environmental matrices. Environ. Sci. Technol. 52, 7409–7417 (2018).CAS 
Article 

Google Scholar 
54.Lusher, A. L., Bråte, I. L. N., Munno, K., Hurley, R. R. & Welden, N. A. Is it or isn’t it: the importance of visual classification in microplastic characterization. Appl. Spectrosc. 74, 1139–1153 (2020).CAS 
Article 

Google Scholar 
55.Primpke, S. et al. Reference database design for the automated analysis of microplastic samples based on Fourier transform infrared (FTIR) spectroscopy. Anal. Bioanal. Chem. 410, 5131–5141 (2018).CAS 
Article 

Google Scholar 
56.Gago, J., Galgani, F., Maes, T. & Thompson, R. Microplastics in seawater: recommendations from the Marine Strategy Framework Directive implementation process. Front. Mar. Sci. https://doi.org/10.3389/fmars.2016.00219 (2016).57.Corcoran, P. L., Belontz, S. L., Ryan, K. & Walzak, M. J. Factors controlling the distribution of microplastic particles in benthic sediment of the Thames River, Canada. Environ. Sci. Technol. 54, 818–825 (2020).CAS 
Article 

Google Scholar  More

Global value of controlling eutrophicationThe substantial emissions from lakes and reservoirs and the potential for increased emissions suggest that there is considerable value in improving water quality in lakes and reservoirs and in preventing further deterioration. We calculated the global climate damages from CH4 emissions and the avoided damages from preventing increased emissions from 2015 to 2050 using well-accepted integrated assessment models (IAMs) (see “Methods”). Because GHGs rapidly become well mixed in the atmosphere, the global social costs of GHG emissions do not depend on where they are emitted. Because GHGs can persist for many years in the atmosphere, the effect of emissions of today will be felt for many years in the future, which means that the rate used to discount future economic damages to the present exerts a strong influence on the social cost of GHG (SC-GHG) estimates. Following the U.S. Government Interagency Working Group (IWG), we report all results using three discount rates: 2.5%, 3%, and 5% yr−1.The estimated present value of the global climate change costs of CH4 emissions from lakes and reservoirs for 2015–2050 range from $7.5 to 81 trillion (2015$; top half of Table 1). Low-end estimates assume a high discount rate (5% yr−1), low current emissions (4.8 Pg CO2-eq yr−1), and no emission growth. High-end estimates assume a low discount rate (2.5% yr−1), high current emissions (8.4 Pg CO2-eq yr−1), and high growth in emissions from lakes (100%). It will not be possible to avoid all emissions from lakes and reservoirs, but with concerted effort it may be possible to prevent increased emissions. The present value of avoided damages from holding emissions constant at current levels rather than increasing by 20–100% by 2050 from increasing eutrophication is $0.66–24 trillion (2015$).Although it has been noted that it might result in underestimation, especially when assuming a high discount rate15, an alternative approach to estimating the climate change damages from non-CO2 GHGs involves first converting the emissions to CO2-equivalents (CO2-eq)16 and then multiplying these by the social cost of carbon dioxide (SC-CO2)15. This approach is less accurate than direct application of the social cost of CH4 (SC-CH4)15, but it has been frequently used in previous studies. To facilitate comparison to other estimates of climate damages in the literature, we also used the CO2-e × SC-CO2 approach with otherwise equivalent assumptions to value eutrophication emissions. Results using this approach are reported in the bottom half of Table 1. The cost of CH4 emissions from lakes and reservoirs from 2015 to 2050 is estimated to be $5.4–95 trillion (2015$), and the associated avoided damages from keeping emissions constant are $0.46–27 trillion (2015$).These estimates consider only the cost of CH4 emissions, but lakes and reservoirs also emit CO2 and N2O. Adding current CO2 and N2O emission estimates10, the SC-GHG emissions increases by 27–51% above those for CH4 alone. Although mounting evidence suggests poor water quality also influences emissions of CO2 and N2O, global analyses of future scenarios for altered emissions of CO2 and N2O from lakes have not yet been published, so we do not monetize these damages. Nevertheless, even our partial estimates suggest that reducing eutrophication is an important means of avoiding climate change damages with a large benefit when measured in monetary terms.Comparison to other economic damages from water pollutionHow do these estimated global climate damages from eutrophication compare to the local and regional benefits of water pollution control typically included in assessments of the benefits and costs of water pollution policies? To help put our results in context, we consider the case of Lake Erie, where eutrophication and associated harmful algal blooms (HABs), primarily due to excess P from agricultural sources, have caused considerable economic damage since the mid-1990s7. Local values of eutrophication abatement vary among lakes, but Lake Erie is a salient example because reliable estimates of local value are available, and Lake Erie’s GHG emissions were included in the global emission analysis9,10 that we used to compute our global estimates presented in Table 1. Recent work using a stated preference survey of Ohio residents estimates that a 40% reduction in total P loading to the western Lake Erie basin from the Maumee River watershed would lead to a $4.0–6.0 million annual welfare gain to Ohio recreational anglers17,18. Assuming constant annual benefits from 2015 to 2050 and using a 3% yr−1 discount rate, this amounts to a present value of $0.087–0.12 billion in total recreational fishing benefits.Table 1 Present value (PV) of global social costs of CH4 emissions from lakes and reservoirs, 2015–2050 (billion 2015 US$).Full size tableApplying our methods to this case, a 40% reduction in total P loading to Lake Erie would yield a 0.079 Tg yr−1 reduction in CH4 emissions (2.7 Tg CO2-eq yr−1). If the P-loading reduction began in 2015 and was maintained through 2050, we estimate that the resulting water quality improvement would generate present value economic benefits (avoided climate damages) of $3.1 billion using the SC-CH4 ($3.3 billion using CO2-e × SC-CO2) and a 3% yr−1 discount rate (Table 2). Thus, the global climate benefits of achieving the targeted 40% reduction in P loading are well over an order of magnitude larger than the estimated recreational benefits to Ohio anglers (Fig. 1).Table 2 Present value (PV) of avoided global social costs of CH4 emissions, 2015–2050 (billion 2015 US$), from a 40% reduction in total P loading in the western Lake Erie basina.Full size tableFig. 1: Comparison of the recreational vs. climate implications of eutrophication.A The welfare gain, 2015–2050, from a 40% reduction in phosphorus (P) loading to western Lake Erie reducing the frequency and extent of harmful algal blooms (HABs). The range of economic impact on recreational angling was estimated from the annual welfare gain17 assuming constant annual benefits and a 3% yr−1 discount rate. The welfare gain from this same total P loading to Lake Erie was estimated from the corresponding reduction in CH4 emissions (and CO2-equivalent emissions) through 2050, using estimates and methods reported in Table 2. B The welfare cost of seasonal Lake Erie HABs sufficient to close beaches, 2015–2050. Benefit transfer work20 estimates the 95% confidence interval of daily recreational losses from the closure of all 67 Lake Erie beaches in Ohio and Michigan. We aggregate to seasonal (115 day)39 HABs occurring annually, 2015–2050, using a 3% yr−1 discount rate. Methane cost estimates are derived from methane emissions under nutrient concentrations that would lead to closure of all of these beaches due to high chlorophyll from HABs as well as from chlorophyll levels that would lead to moderate risk of adverse health effects from beach use.Full size imagePublished estimates suggest that the 40% reduction in total P loading to Lake Erie that we model here could be achieved with a fertilizer tax or a tax-and-rebate policy with rebates funding agricultural best management practices at an annual cost to taxpayers of about $16–17 million19. Note that these cost estimates are conservative, as they do not include yield losses or other agricultural compliance costs19. These annual costs would exceed the estimated annual recreational fishing benefits of the policy goal18 but are still smaller than the climate benefits.Economists have also used benefit transfer techniques to extrapolate from individual estimates of the value of water quality changes for a specific location to estimates for an entire region. For example, recent work20 using a function transfer approach estimates that the closure of all 67 Lake Erie beaches in Ohio and Michigan due to a large HAB in Lake Erie would generate daily recreational losses of $2.39 million (95% confidence interval $1.81–3.11 million). Assuming an extreme case that the HAB season lasts continuously for 115 days20, this implies an annual welfare loss of about $280 million. If a severe HAB that closed all 67 Lake Erie beaches in the two states occurs annually from 2015 to 2050 and annual damages are constant, the present value of total damages, derived from the definition of the present value of a constant stream of benefits, using a 3% yr−1 discount rate, would be about $6.1 billion using the central estimate of the cost of beach closure20, or a range of $4.4–7.7 billion, using their 95% confidence interval20.The CH4 emissions from a HAB event in Lake Erie large enough to close all 67 beaches in Ohio and Michigan would depend on the severity of the triggering water quality impairment. We use two approaches to make a comparable estimate of CH4 emission damages. First, if the chlorophyll a concentration exceeds 30 ppb, the risk of Cyanobacteria blooms is 80–100%, gauged by the risk of Cyanobacteria biomass exceeding 50%21. This level would exceed statutory thresholds that trigger beach closures or health advisories and would yield an emission increase of 1.0 Tg CH4 yr−1 (34 Tg CO2-eq yr−1). These emissions would create a present value of damages of $39 billion using the SC-CH4 ($42 billion using CO2-e × SC-CO2) at a 3% yr−1 discount rate (Table 3), roughly seven times larger than the estimated recreational damages from a HAB severe enough to close all Lake Erie beaches in Michigan and Ohio for 35 years.Table 3 Present value (PV) of global social costs of CH4 emissions, 2015–2050 (billion 2015 US$), from a harmful algal bloom sufficient to close all MI and OH beaches on Lake Erie.Full size tableAs a second approach to making this comparison, we use the World Health Organization guideline for chlorophyll a concentration yielding a moderate probability of adverse health effects in recreational waters (50 ppb)22. Because the assumed triggering concentration for beach closures is higher, both the estimated emissions associated with the closure events (1.7 Tg CH4 yr−1 or 59 Tg CO2-eq yr−1) and the economic damages using a 3% yr−1 discount rate ($69 and $73 billion) are higher (Table 3). With this approach, the global climate costs of HABs severe enough to close all MI and OH beaches on Lake Erie from 2015 to 2050 are an order of magnitude larger than the estimated recreational damages from beach closures (Fig. 1).We cannot say how our CH4 damage estimates would compare with a full estimate of other damages from Lake Erie eutrophication. The literature demonstrates that important water quality benefits are difficult to value2. A single-season HAB similar to the 2014 event that resulted in the issuance of a do not drink/do not boil order for the public water system in the City of Toledo created damages of about $1.3 billion, including impacts on property values, water treatment costs, and tourism23. Estimates of damages to fishing activity at Lake Erie’s Canadian coast are also substantial24. An earlier study estimates damages from eutrophication of all U.S. rivers and lakes25, omitting the climate damage estimates we calculate here; an assessment of the methods used to obtain these estimates is outside the scope of our paper. Notably, recent work links HABs in Gull Lake, Michigan (not far from Lake Erie) with increased likelihood of low birth weight and shorter gestation among infants born to exposed mothers26.Given that the full gamut of potential damages is difficult to monetize, a comprehensive estimate of the non-climate damages from eutrophication and HABs—especially if human health impacts are significant—could exceed our damage estimates for CH4 emissions. However, our estimates of the global CH4 emission damages from eutrophication in Lake Erie exceed all published estimates of other damages, to the extent that we can compare them. Smaller lakes than Erie may show even greater differences between global and local values of eutrophication because, on average, people have greater willingness to pay for recreation on large lakes27, and CH4 emissions per unit area do not vary with lake size10. These results suggest that global climate impacts are a substantial omission from benefit–cost assessments of policies targeting eutrophication, in Lake Erie and elsewhere.Eutrophication is a local and global problemDegraded water quality is often considered a local or regional problem. We show that water quality has important implications for global climate, through emissions of CH4 and other GHGs. These emissions are likely to increase substantially unless action is taken to prevent further eutrophication. The damage from eutrophication-related GHG emissions is likely to be in trillions of dollars, and appears to be far larger than other monetized damages from poor water quality that economists have so far been able to quantify, especially where pollution does not generate severe health damages. Our analysis shows that local water quality protection has global economic implications, and that more effort devoted to understanding the consequences of changes in water quality and valuing the benefits of sustaining or improving water quality is warranted. More

1.Drought in Central‐Northern Europe (EDO, 2018); https://edo.jrc.ec.europa.eu/documents/news/EDODroughtNews201809_Central_North_Europe.pdf2.Drought in Europe (EDO, 2019); https://edo.jrc.ec.europa.eu/documents/news/EDODroughtNews201908_Europe.pdf3.Kovats, R. S. et al. in Climate Change 2014: Impacts, Adaptation, and Vulnerability (eds Fields, C. B. et al.) 1267–1326 (Cambridge Univ. Press, 2014).4.Spinoni, J., Naumann, G. & Vogt, J. V. Pan-European seasonal trends and recent changes of drought frequency and severity. Glob. Planet. Change 148, 113–130 (2017).Article 

Google Scholar 
5.Vörösmarty, C. J., Green, P., Salisbury, J. & Lammers, R. B. Global water resources: vulnerability from climate change and population growth. Science 289, 284–288 (2000).Article 

Google Scholar 
6.Döll, P., Fiedler, K. & Zhang, J. Global-scale analysis of river flow alterations due to water withdrawals and reservoirs. Hydrol. Earth Syst. Sci. 13, 2413 (2009).Article 

Google Scholar 
7.Wada, Y., Van Beek, L. P., Wanders, N. & Bierkens, M. F. Human water consumption intensifies hydrological drought worldwide. Environ. Res. Lett. 8, 034036 (2013).Article 

Google Scholar 
8.Tijdeman, E., Hannaford, J. & Stahl, K. Human influences on streamflow drought characteristics in England and Wales. Hydrol. Earth Syst. Sci. 22, 1051–1064 (2018).Article 

Google Scholar 
9.Beniston, M. et al. Future extreme events in European climate: an exploration of regional climate model projections. Climatic Change 81, 71–95 (2007).Article 

Google Scholar 
10.Nikulin, G., Kjellstrom, E., Hansson, U. L. F., Strandberg, G. & Ullerstig, A. Evaluation and future projections of temperature, precipitation and wind extremes over Europe in an ensemble of regional climate simulations. Tellus A 63, 41–55 (2011).Article 

Google Scholar 
11.Forzieri, G. et al. Ensemble projections of future streamflow droughts in Europe. Hydrol. Earth Syst. Sci. 18, 85–108 (2014).Article 

Google Scholar 
12.Samaniego, L. et al. Anthropogenic warming exacerbates European soil moisture droughts. Nat. Clim. Change 8, 421 (2018).Article 

Google Scholar 
13.Marx, A. et al. Climate change alters low flows in Europe under global warming of 1.5, 2, and 3 °C. Hydrol. Earth Syst. Sci. 22, 1017–1032 (2018).Article 

Google Scholar 
14.Stahl, K. et al. Impacts of European drought events: insights from an international database of text-based reports. Nat. Hazards Earth Syst. Sci. 16, 801–819 (2016).Article 

Google Scholar 
15.Mapping the Impacts of Natural hazards and Technological Accidents in Europe: An Overview of the Last Decade Technical Report No. 13/2010 (EEA, 2011); http://op.europa.eu/en/publication-detail/-/publication/4f5878ba-0947-4fb6-964b-8818cfda3de7/language-en16.Schär, C. et al. The role of increasing temperature variability in European summer heatwaves. Nature 427, 332–336 (2004).Article 
CAS 

Google Scholar 
17.Gil, M., Garrido, A. & Hernández-Mora, N. Direct and indirect economic impacts of drought in the agri-food sector in the Ebro River basin (Spain). Nat. Hazards Earth Syst. Sci. 13, 2679–2694 (2013).Article 

Google Scholar 
18.García-León, D., Standardi, G. & Staccione, A. An integrated approach for the estimation of agricultural drought costs. Land Use Policy 100, 104923 (2021).Article 

Google Scholar 
19.Byers, E. A., Coxon, G., Freer, J. & Hall, J. W. Drought and climate change impacts on cooling water shortages and electricity prices in Great Britain. Nat. Commun. 11, 2239 (2020).CAS 
Article 

Google Scholar 
20.Salmoral, G., Rey, D., Rudd, A., Margon, Pde & Holman, I. A probabilistic risk assessment of the national economic impacts of regulatory drought management on irrigated agriculture. Earths Future 7, 178–196 (2019).Article 

Google Scholar 
21.Naumann, G., Spinoni, J., Vogt, J. V. & Barbosa, P. Assessment of drought damages and their uncertainties in Europe. Environ. Res. Lett. 10, 124013 (2015).Article 

Google Scholar 
22.Stagge, J. H., Kohn, I., Tallaksen, L. M. & Stahl, K. Modeling drought impact occurrence based on meteorological drought indices in Europe. J. Hydrol. 530, 37–50 (2015).Article 

Google Scholar 
23.Blauhut, V. et al. Estimating drought risk across Europe from reported drought impacts, drought indices, and vulnerability factors. Hydrol. Earth Syst. Sci. 20, 2779–2800 (2016).Article 

Google Scholar 
24.Freire-González, J., Decker, C. & Hall, J. W. The economic impacts of droughts: a framework for analysis. Ecol. Econ. 132, 196–204 (2017).Article 

Google Scholar 
25.The 2015 Ageing Report: Underlying Assumptions and Projection Methodologies (European Commission, 2014).26.Berg, A. et al. Land–atmosphere feedbacks amplify aridity increase over land under global warming. Nat. Clim. Change 6, 869–874 (2016).Article 

Google Scholar 
27.Dosio, A. & Fischer, E. M. Will half a degree make a difference? Robust projections of indices of mean and extreme climate in Europe under 1.5 °C, 2 °C, and 3 °C global warming. Geophys. Res. Lett. 45, 935–944 (2018).Article 

Google Scholar 
28.Jacob, D. et al. Climate impacts in Europe under +1.5 °C global warming. Earths Future 6, 264–285 (2018).Article 

Google Scholar 
29.Alfieri, L., Dottori, F., Betts, R., Salamon, P. & Feyen, L. Multi-model projections of river flood risk in Europe under global warming. Climate 6, 6 (2018).Article 

Google Scholar 
30.Vousdoukas, M. I. et al. Climatic and socioeconomic controls of future coastal flood risk in Europe. Nat. Clim. Change 8, 776–780 (2018).Article 

Google Scholar 
31.Estrela, T. & Vargas, E. Drought management plans in the European Union. The case of Spain. Water Resour. Manag. 26, 1537–1553 (2012).Article 

Google Scholar 
32.Dellink, R., Chateau, J., Lanzi, E. & Magné, B. Long-term economic growth projections in the shared socioeconomic pathways. Glob. Environ. Change 42, 200–214 (2017).Article 

Google Scholar 
33.Christensen, P., Gillingham, K. & Nordhaus, W. Uncertainty in forecasts of long-run economic growth. Proc. Natl Acad. Sci. USA 115, 5409–5414 (2018).CAS 
Article 

Google Scholar 
34.Global Assessment Report on Disaster Risk Reduction 2019 (United Nations Office for Disaster Risk Reduction, 2019).35.Forzieri, G. et al. Escalating impacts of climate extremes on critical infrastructures in Europe. Glob. Environ. Change 48, 97–107 (2018).Article 

Google Scholar 
36.Erfurt, M., Glaser, R. & Blauhut, V. Changing impacts and societal responses to drought in southwestern Germany since 1800. Reg. Environ. Change 19, 2311–2323 (2019).37.Swiss Re The hidden risks of climate change: an increase in property damage from soil subsidence in Europe. PreventionWeb https://www.preventionweb.net/publications/view/20623 (2011).38.Formetta, G. & Feyen, L. Empirical evidence of declining global vulnerability to climate-related hazards. Glob. Environ. Change 57, 101920 (2019).Article 

Google Scholar 
39.Zhang, H., Li, Y. & Zhu, J.-K. Developing naturally stress-resistant crops for a sustainable agriculture. Nat. Plants 4, 989–996 (2018).Article 

Google Scholar 
40.Lohrmann, A., Farfan, J., Caldera, U., Lohrmann, C. & Breyer, C. Global scenarios for significant water use reduction in thermal power plants based on cooling water demand estimation using satellite imagery. Nat. Energy 4, 1040–1048 (2019).Article 

Google Scholar 
41.Hallegatte, S., Przyluski, V. & Vogt-Schilb, A. Building world narratives for climate change impact, adaptation and vulnerability analyses. Nat. Clim. Change 1, 151–155 (2011).Article 

Google Scholar 
42.Di Baldassarre, G. et al. Water shortages worsened by reservoir effects. Nat. Sustain. 1, 617–622 (2018).Article 

Google Scholar 
43.Guerreiro, S. B., Dawson, R. J., Kilsby, C., Lewis, E. & Ford, A. Future heat-waves, droughts and floods in 571 European cities. Environ. Res. Lett. 13, 034009 (2018).Article 

Google Scholar 
44.Vetter, T. et al. Evaluation of sources of uncertainty in projected hydrological changes under climate change in 12 large-scale river basins. Climatic Change 141, 419–433 (2017).CAS 
Article 

Google Scholar 
45.Hattermann, F. F. et al. Sources of uncertainty in hydrological climate impact assessment: a cross-scale study. Environ. Res. Lett. 13, 015006 (2018).Article 

Google Scholar 
46.Hattermann, F. F. et al. Cross-scale intercomparison of climate change impacts simulated by regional and global hydrological models in eleven large river basins. Climatic Change 141, 561–576 (2017).Article 

Google Scholar 
47.Addressing the Challenge of Water Scarcity and Droughts in the European Union (European Commission, 2007); https://www.eea.europa.eu/policy-documents/addressing-the-challenge-of-water48.Smith, A. B. U.S. Billion-Dollar Weather and Climate Disasters, 1980–Present (NCEI Accession 0209268) (NOAA, 2020); https://doi.org/10.25921/STKW-7W7349.Martin-Ortega, J., González-Eguino, M. & Markandya, A. The costs of drought: the 2007/2008 case of Barcelona. Water Policy 14, 539–560 (2012).Article 

Google Scholar 
50.Zampieri, M. et al. Climate resilience of the top ten wheat producers in the Mediterranean and the Middle East. Reg. Environ. Change 20, 41 (2020).Article 

Google Scholar 
51.Vliet, M. T. H., van, Vögele, S. & Rübbelke, D. Water constraints on European power supply under climate change: impacts on electricity prices. Environ. Res. Lett. 8, 035010 (2013).Article 

Google Scholar 
52.Lehner, B., Czisch, G. & Vassolo, S. The impact of global change on the hydropower potential of Europe: a model-based analysis. Energy Policy 33, 839–855 (2005).Article 

Google Scholar 
53.Jenkins, K. Indirect economic losses of drought under future projections of climate change: a case study for Spain. Nat. Hazards 69, 1967–1986 (2013).Article 

Google Scholar 
54.Gall, M., Borden, K. A. & Cutter, S. L. When do losses count? Six fallacies of natural hazards loss data. Bull. Am. Meteorol. Soc. 90, 799–810 (2009).Article 

Google Scholar 
55.Xu, C. et al. Increasing impacts of extreme droughts on vegetation productivity under climate change. Nat. Clim. Change 9, 948–953 (2019).CAS 
Article 

Google Scholar 
56.Choat, B. et al. Triggers of tree mortality under drought. Nature https://www.nature.com/articles/s41586-018-0240-x (2018).57.Seidl, R. et al. Invasive alien pests threaten the carbon stored in Europe’s forests. Nat. Commun. 9, 1626 (2018).Article 
CAS 

Google Scholar 
58.Sutanto, S. J., Vitolo, C., Di Napoli, C., D’Andrea, M. & Van Lanen, H. A. J. Heatwaves, droughts, and fires: exploring compound and cascading dry hazards at the pan-European scale. Environ. Int. 134, 105276 (2020).Article 

Google Scholar 
59.de Ruiter, M. C. et al. Why we can no longer ignore consecutive disasters. Earths Future 8, e2019EF001425 (2020).Article 

Google Scholar 
60.Ford, T. W. & Labosier, C. F. Meteorological conditions associated with the onset of flash drought in the eastern United States. Agric. For. Meteorol. 247, 414–423 (2017).Article 

Google Scholar 
61.Yuan, X., Wang, L. & Wood, E. F. Anthropogenic intensification of southern African flash droughts as exemplified by the 2015/16 season. Bull. Am. Meteorol. Soc. 99, S86–S90 (2018).Article 

Google Scholar 
62.Yuan, X., Ma, Z., Pan, M. & Shi, C. Microwave remote sensing of short-term droughts during crop growing seasons. Geophys. Res. Lett. 42, 4394–4401 (2015).Article 

Google Scholar 
63.Nguyen, H. et al. Using the evaporative stress index to monitor flash drought in Australia. Environ. Res. Lett. 14, 064016 (2019).Article 

Google Scholar 
64.Pendergrass, A. G. et al. Flash droughts present a new challenge for subseasonal-to-seasonal prediction. Nat. Clim. Change 10, 191–199 (2020).Article 

Google Scholar 
65.Hagenlocher, M. et al. Drought vulnerability and risk assessments: state of the art, persistent gaps, and research agenda. Environ. Res. Lett. 14, 083002 (2019).Article 

Google Scholar 
66.Jacobs-Crisioni, C. et al. The LUISA Territorial Reference Scenario 2017: A Technical Description (Publications Office of the European Union, 2017); https://ec.europa.eu/jrc/en/publication/luisa-territorial-reference-scenario-201767.Capros, P. et al. GEM-E3 Model Documentation (Publications Office of the European Union, 2013); https://publications.jrc.ec.europa.eu/repository/handle/111111111/3236668.Keramidas, K., Kitous, A., Després, J. & Schmitz, A. POLES-JRC Model Documentation (Publications Office of the European Union, 2017); https://publications.jrc.ec.europa.eu/repository/handle/JRC11375769.Feyen, L. & Dankers, R. Impact of global warming on streamflow drought in Europe. J. Geophys. Res. Atmos. 114, D17116 (2009).Article 

Google Scholar 
70.Tallaksen, L. M. & Van Lanen, H. A. Hydrological Drought: Processes and Estimation Methods for Streamflow and Groundwater Vol. 48 (Elsevier, 2004).71.Lehner, B., Döll, P., Alcamo, J., Henrichs, T. & Kaspar, F. Estimating the impact of global change on flood and drought risks in Europe: a continental, integrated analysis. Climatic Change 75, 273–299 (2006).Article 

Google Scholar 
72.Roudier, P. et al. Projections of future floods and hydrological droughts in Europe under a +2 °C global warming. Climatic Change 135, 341–355 (2016).Article 

Google Scholar 
73.Knijff, J. M. V. D., Younis, J. & Roo, A. P. J. D. LISFLOOD: a GIS‐based distributed model for river basin scale water balance and flood simulation. Int. J. Geogr. Inf. Sci. 24, 189–212 (2010).Article 

Google Scholar 
74.Salamon, P. et al. EFAS Upgrade for the Extended Model Domain (Publications Office of the European Union, 2019); https://publications.jrc.ec.europa.eu/repository/handle/111111111/5558775.Jacob, D. et al. EURO-CORDEX: new high-resolution climate change projections for European impact research. Reg. Environ. Change 14, 563–578 (2014).Article 

Google Scholar 
76.Knutti, R. et al. Meeting Report. In IPCC Expert Meeting on Assessing and Combining Multi Model Climate Projections (eds Stocker, T. F. et al.) (IPCC, 2010).77.Shepherd, T. G. Storyline approach to the construction of regional climate change information. Proc. R. Soc. Lond. A 475, 20190013 (2019).
Google Scholar 
78.Mentaschi, L. et al. Independence of future changes of river runoff in Europe from the pathway to global warming. Climate 8, 22 (2020).Article 

Google Scholar 
79.Mentaschi, L. et al. The transformed-stationary approach: a generic and simplified methodology for non-stationary extreme value analysis. Hydrol. Earth Syst. Sci. 20, 3527–3547 (2016).Article 

Google Scholar 
80.Forzieri, G. et al. Resilience of Large Investments and Critical Infrastructures in Europe to Climate Change (Publications Office of the European Union, 2015); https://publications.jrc.ec.europa.eu/repository/handle/111111111/3889481.Batista e Silva, F. et al. HARCI-EU, a harmonized gridded dataset of critical infrastructures in Europe for large-scale risk assessments. Sci. Data 6, 126 (2019).Article 

Google Scholar 
82.Doornkamp, J. C. Clay shrinkage induced subsidence. Geogr. J. 159, 196–202 (1993).Article 

Google Scholar 
83.Boivin, P., Garnier, P. & Tessier, D. Relationship between clay content, clay type, and shrinkage properties of soil samples. Soil Sci. Soc. Am. J. 68, 1145–1153 (2004).CAS 
Article 

Google Scholar 
84.Hiederer, R. Mapping Soil Properties for Europe—Spatial Representation of Soil Database Attributes (Publications Office of the European Union, 2013); https://publications.jrc.ec.europa.eu/repository/handle/111111111/2917085.Crilly, M. Analysis of a database of subsidence damage. Struct. Surv. 19, 7–15 (2001).Article 

Google Scholar 
86.Corti, T., Wüest, M., Bresch, D. & Seneviratne, S. I. Drought-induced building damages from simulations at regional scale. Nat. Hazards Earth Syst. Sci. 11, 3335–3342 (2011).Article 

Google Scholar 
87.Batista e Silva, F., Lavalle, C. & Koomen, E. A procedure to obtain a refined European land use/cover map. J. Land Use Sci. 8, 255–283 (2013).Article 

Google Scholar 
88.Florczyk, A. et al. GHSL Data Package 2019 (Publications Office of the European Union, 2019); https://publications.jrc.ec.europa.eu/repository/handle/111111111/5655289.Kron, W., Steuer, M., Löw, P. & Wirtz, A. How to deal properly with a natural catastrophe database—analysis of flood losses. Nat. Hazards Earth Syst. Sci. 12, 535–550 (2012).90.Felbermayr, G. & Gröschl, J. Naturally negative: the growth effects of natural disasters. J. Dev. Econ. 111, 92–106 (2014).Article 

Google Scholar  More

Quantitative variablesThe intervention described here relies on the IT platform “SmartH2O” for the collection and visualization of smart meter data, the provision of consumption feedback to the user, the delivery of water-saving recommendations, and the engagement of the consumer through a gamification program20,21,34,35. We embedded a gamification mechanism in the digital platform to maximize user retention and stimulate the exploration and sharing of content and the setting and achievement of personal saving goals. Via the gamification mechanisms, users could collect reward points for different actions performed in the digital platform or the achievement of water-saving targets. Reward points consisted of virtual points that the users could redeem for physical rewards. The design of the SmartH2O digital platform and the behavioral change stimuli that have been introduced in the Valencia case-control study (e.g., web and mobile app, different reward schemes), along with their individual elements and the corresponding illustrative screenshots of the platform are provided in the Supplementary Information, consistently with the information published in a previous study21. Other platforms similar to SmartH2O or approaches for water conservation based on digital technologies are reported in the literature, including, e.g., real-time water consumption feedback on in-home displays, interactive dashboards, and games36,37. Yet, to the author’s knowledge, SmartH2O is the first platform of its kind whose effect is rigorously assessed in the medium term and long term.Household-scale water consumption data and smart meter sampling frequencyWater consumption readings measured at the household scale constitute the main quantitative variable of interest used in this observational study to identify behavior changes. The SmartH2O digital platform relies on water consumption information stored in a central database and enables data communication from the water utility to the water consumers (see Supplementary Fig. 1 for its software architecture). Water consumption data are collected by smart meters installed at the household premises, according to a schedule that considers the maximum available frequency of data sampling at each installation (hourly or daily). The consumption data are anonymized by the utility company, filtered, and transferred to the central database of the SmartH2O platform. The content of the central database is published to the user via a web portal and a mobile application, which are the entry points of all users’ interactions with the platform.Besides the time series of water consumption, we also stored the sampling frequency allowed by each household-scale smart meter. Two types of sampling frequencies were available in the considered population, depending on the installed smart meter hardware: hourly or daily.Digital user engagement variablesThe central database of the SmartH2O platform comprises content for improving user awareness, such as water-saving recommendations, and for implementing the gamification program, such as the description of virtual and physical rewards. The interaction of the users with the platform and the overall user experience features several functionalities, including user login, water consumption and smart meter-based feedback visualization, conservation goal settings, and different gamified water conservation awareness actions (see also Supplementary Notes1). We monitored the activity of each user in the SmartH2O platform for the entire duration of the treatment period and gathered quantitative data on these four digital user engagement variables:

(i)

Login count, defined as the total number of logins executed by each user.

(ii)

Non-rewarded action count, defined as the total number of actions performed by each user, with no reward points associated.

(iii)

Rewarded action count, defined as the total number of actions performed by each user, with associated reward points upon their completion.

(iv)

Cumulative reward importance, defined as the total amount of points achieved by each user by completing the rewarded actions. It accounts for the total amount of points, badges, and rewards achieved by an individual user in the SmartH2O platform.

Each user profile in the SmartH2O platform was associated with a unique smart meter ID, which allowed linking the user activity in the platform with the household water consumption data. User confidentiality was maintained throughout the full study as data were anonymized by the water utility managing the water meters and the central database.Population and study sizeOur observational study was conducted in the city of Valencia, Spain. With a population of 794,288 inhabitants, as reported in 2019 by the Spanish National Institute of Statistics (Institudo Nacional de Estadística)38, Valencia is the third-largest city in Spain. The water utility of Valencia (Global Omnium–EMIVASA) has installed more than 425,000 smart water meters since the early developments in 2006 to monitor the water consumption of nearly all the population39 (the last official census data, recorded in 2011, report 419,994 households in total in Valencia40). The total population considered in this study after application of the exclusion criteria described in the next section included 334 individual households, each equipped with a water meter.The architecture of the smart metering infrastructure deployed in Valencia has been designed in order to be vendor-independent, so it allows for different smart metering solutions to be integrated39. While this is clearly an advantage for procurement, the diversity of hardware has an impact on data sampling and only one of the available technologies supports hourly data collection, which is a preferred requirement for water consumption data quality assessment and provision of sub-daily water consumption information to households in our case-control study. The number of hourly reading meters in Valencia amounts to 168,172 as of July 12th, 2020. EMIVASA also offered its customers access to a web platform where bills and invoices could be managed and also information about the current (daily and monthly) water consumption data was made available.During our observational study, we integrated the digital SmartH2O platform20,21 in the EMIVASA portal. We invited users who already had an account in the platform and a compatible meter reading frequency to voluntarily join our observational study and sign up to the SmartH2O platform. The recruitment campaign was performed using different media channels, namely, newspaper articles on consumer magazines, radio programs, banners on the digital and printed invoices sent to EMIVASA customers, and also a Facebook campaign targeting the Valencia area. At the end of the recruitment campaign, we received 525 applications out of which we obtained a treatment group composed of 223 households after application of the inclusion/exclusion criteria. Out of the households who did not apply to join the case-control study during the recruitment phase, 111 households agreed to be monitored as part of the self-selected control group to be considered as a benchmark group not subject to treatment, after active recruitment via phone by the EMIVASA call center (client service management). Households in the control group had only access to their water consumption data through the already existing platform, which did not offer any type of smart meter-based consumption feedback, behavioral stimuli, and/or gamification elements.Informed consent was obtained from the households monitored in this study. Moreover, the water utility (Global Omnium–EMIVASA) supervised and approved the collection, usage, and processing of the anonymized quantitative variables above described in compliance with the EU General Data Protection Regulation 2016/679 and the pre-existing Spanish law 15/1999 LOPD of 1999 (the SmartH2O study started before the adoption of the GDPR in 2016).Baseline and observation periodsThe treatment period of the case-control study lasted 8.5 months, from June 2016 to February 2017. We also continuously collected anonymized water consumption data for the study population from June 2016 to February 2019 both to conduct the longitudinal study presented in this paper and evaluate water consumption changes over time in comparison with a pre-treatment baseline (June 1st, 2015– April 30th, 2016), as well as to compare water consumption changes in the treatment and control groups. Consistently with the months included in the treatment period (short-term behavior change), we identify the observation period June 1st, 2017–February 2nd, 2018 for medium-term behavior change assessment, and the observation period June 1st, 2018–February 2nd, 2019 for long-term behavior change.Exclusion criteriaThe population considered for analysis of water consumption changes in this observational study was obtained by sequential application of the following exclusion criteria.

1.

Exclusion of empty households. First, we excluded the households with no data in the baseline and treatment period. We classified in this category also the households with a cumulative water consumption lower than 1.5 m3 over the whole baseline and treatment period (which together last nearly 20 months). This threshold value was identified as a conservative choice after consultation with the local water utility and comparison with the average values of water consumption in the entire population (slightly above 0.21 m3/day) and the European average water consumption, which amounts to 128 liters per inhabitant per day (0.128 m3/day)41. A household in the considered population would use ~1.5 m3 in one week (0.21 m3/day × 7 days). While lower values than the average consumption are observed in those days in which the inhabitants spend little time at home, a cumulative consumption of 1.5 m3 over the course of more than 1 year can indicate that the house is generally empty (and possibly the observed water consumption is due to leaks).

2.

Exclusion of households with insufficient data length. We removed the households with water consumption readings for less than 1000 h (approximately 6 weeks). This step guarantees a minimum representation of weekend/weekday water demand variation for more than 1 month (please note that the total duration of the treatment period is 8.5 months).

3.

Exclusion of partially empty households. We excluded the households with more than 90% water consumption readings equal to zero in the baseline or observation period or completely lacking data for one of these two periods. We considered these households to be empty or equipped with faulty meters at least during one of the two short-term periods of interest. The above value threshold of 90% was identified with a trial-and-error procedure and expert-based data analysis that balance the rate of exclusion with the size of the remaining dataset.

4.

Exclusion of households lacking day-of-week representation. We excluded the households with available observations for less than 7 unique day types, to guarantee a minimum representation of water consumption routines that depend on the day of the week. For those households with smart meters recording water consumption with hourly sampling frequency, we removed days with more than 4 h of gaps from the smart meter time series (anomalous meter data logging).

5.

Exclusion of households with anomalous high water consumption. We considered hourly water consumption readings larger than 1 m3 as outliers (we thus removed these hourly readings) and we removed the households with a daily average water consumption larger than 1 m3 in at least one phase of the longitudinal study. High values of water consumption can be observed for specific days (e.g., when customers use water for outdoor irrigation or filling up a pool), yet average daily water consumption values over the selected threshold are more than three times higher than the European average (equivalent to approximately 0.3 m3/day per household). We did not apply more restrictive thresholds, in order not to bias our analysis and avoid unjustified exclusion of high water consumers.

6.

Exclusion of households with unrealistic short-term consumption change levels. We excluded the households with extreme values of short-term consumption change during the treatment period, which were identified as outliers by Tukey’s fences42. According to Tukey’s fences, a data point xi is considered an outlier if:$$x_i notin [Q_1 – kleft( {Q_3 – Q_1} right),Q_3 + k(Q_3 – Q_1)]$$
(1)

where Q1 is the 25th empirical quartile (i.e., 25% of the data is lower than this point) and Q3 is the 75th empirical quartile (i.e., 75% of the data is lower than this point), and k = 1.5. Tukey’s fences with k = 1.5 approximate the 99.7% confidence interval defined for normal distributions by a distance of three standard deviations from the mean.

7.

Exclusion of households with anomalous conditions in medium-term and long-term. We excluded 51 households that met the above exclusion criteria 1–6 during either the medium-term or long-term observation periods. Water consumption change patterns would be incomplete/anomalous for these households, with at least one missing/anomalous period out of the four periods of interest (i.e., baseline, treatment period, or following observation periods in 2018 and 2019).

With the above exclusion criteria, we obtained the 334 households considered for behavior change analysis in this observational study. More details on the population size after application of each exclusion criteria are reported with a flow diagram in Supplementary Fig. 1043. It is worth noting that only less than 2% of high consumption households have been excluded, while most of the other excluded households had insufficient data or unrealistically low consumption levels. Also, the number of households in the sample considered here differ from those considered in the evaluation of the SmartH2O project44, due to the different temporal length of the two studies and the application of the exclusion criteria on data recorded in different periods (the SmartH2O project only included the baseline and treatment periods).Adopting the same criteria to exclude households from the behavior change analysis only during the summer period (Fig. 1d) resulted in a reduced population of 179 households (101 households in the treatment group and 78 households in the control group), due to limited data availability for the summer period. Similarly, a subset of 198 households in the treatment group was considered for the correlation analysis by logistic regression (Fig. 4), as the excluded 25 households presented incomplete smart meter data or incomplete information on their usage of the digital SmartH2O application.Data analysis and statistical methodsWe performed customer segmentation to analyze heterogeneous long-term behavior change patterns (Fig. 3 and Supplementary Fig. 9). We applied agglomerative hierarchical clustering45 to the patterns of average daily household water consumption during the entire duration of the longitudinal study. Here, a water consumption pattern of a household is a vector that contains four values of average daily water consumption, i.e., one for each period of the observational study, including the baseline (see “Methods” section–Baseline and observation periods). The only variable given as input to the hierarchical clustering algorithm consists of household-scale average water consumption per day for each phase of our observational study, which spans the baseline and the three observation periods in 2017, 2018, and 2019. Complete linkage and correlation distance were considered for hierarchical clustering. Complete linkage calculates the distance between two household clusters as the distance between the farthest pair of household water consumption patterns in the two household clusters, i.e., the maximum distance formulated as follows:46$$dleft( {u,v} right) = {mathrm{max}}left( {{mathrm{dist}}left( {uleft( {x_i} right),vleft( {z_i} right)} right)} right.$$
(2)
where d(u,v) is the distance between clusters u and v, xi are the points belonging to cluster u and zi those belonging to cluster v. Given two vectors of observations xi and xj, which in our study correspond to the water consumption patterns of two households (each with N elements, with N = 4, where each element is the household-scale average water consumption per day for the baseline and three observation periods) and their mean values ((bar x_i) and (bar x_j)) the correlation distance used by hierarchical clustering is calculated as follows:47$${mathrm{dist}}left( {x_i,x_j} right) = 1 – frac{{(x_i – bar x_i) cdot (x_j – bar x_j)}}{{left| {x_i – bar x_i} right|_2left| {x_j – bar x_j} right|_2}}$$
(3)
We considered hierarchical clustering as an appropriate choice because the analysis of the different hierarchical levels allowed the discovery of heterogeneous water consumption behaviors that would be potentially hidden if algorithms requiring a predefined number of clusters were used. We adopted complete linkage clustering to avoid that individual, mutually close households would force pairs of clusters representing different behaviors to merge. Also, we adopted correlation distance as we wanted to identify similarities in water consumption patterns over time, rather than in water consumption volumes.After clustering the households in the treatment group with the above hierarchical clustering, similarly to a previous study18, we analyzed the coefficients of a logistic regression classifier cross-validated with binary tests to identify which candidate factors correlate with the main behavior change patterns that characterize the households in the treatment group (Fig. 4 and Supplementary Table 3). In this study, the input candidate factors consist of five independent variables that comprise the availability of smart meter hourly data frequency and the four digital user engagement variables, i.e., login count, non-rewarded action count, rewarded action count, and cumulative reward importance. First, we balanced the distribution of the households in the treatment group across the behavior change segments considered in the binary tests by Synthetic Minority Over-sampling Technique (SMOTE)48. SMOTE oversamples the minority class to balance the sample distribution of a labeled dataset over the different classes. As we consider binary test where only two behavior change segments (or two groups of behavior change segments) are compared, the majority class represents the behavior change segment (or group of behavior change segments) with the highest number of samples and vice versa for the minority class. According to the SMOTE formulation48, starting from a sample ci,initial, which in this study is the vector of input candidate factors for a household i in the minority class, a new sample ci,new is generated on the line between ci,initial, and one of its k nearest-neighbors cj,initial, with the following formula:$$c_{i,{mathrm{new}}} = c_{i,{mathrm{initial}}} + lambda (c_{j,{mathrm{initial}}} – c_{i,{mathrm{initial}}})$$
(4)
where λ is a random number between 0 and 1, and k = 5 nearest neighbors computed based on Euclidean distance are considered by default48. Among the possible options to perform class balancing, here we adopted a “not majority” strategy to over-sample the minority classes, i.e., we resample all classes but the majority class (which, in our binary problem, is equivalent to resampling the minority class).Second, we trained a logistic regression classifier49 with k-fold cross-validation (k = 5) and evaluated its performance via weighted F1 score. In our binary problem, the logistic regression classifier models the class membership probability P(yi,p = 1) for household i, where yi,p = 1 indicates that the household belongs to behavior change pattern p (else yi,p = 0, according to the following logistic function:$$Pleft( {y_{i,p} = 1} right) = frac{1}{{1 + exp^{ – f(c_i)}}}$$
(5)
where f(ci) is a linear function where the input variables ci are weighted by corresponding coefficients α:$$fleft( {c_i} right) = alpha _0 + alpha _1c_{i,1} + alpha _2c_{i,2} + ldots + alpha _Mc_{i,M} + varepsilon _i$$
(6)
In this study, M = 5, ci,1 is a binary variable representing the availability of smart meter with hourly data frequency, ci,{2,3,4,5} are the four digital user engagement variables defined above, α0 is the intercept of the logistic regression, and εi is random noise. We normalized the variables before logistic regression classification by subtracting the mean and dividing by the standard deviation to rescale them to comparable value ranges. The analysis of their corresponding logistic regression coefficients reveals how these variables discriminate among different clusters of water consumers and, thus, how they are potential determinants of defined water consumption behaviors. The F1 score (FS) is first calculated for each behavior change pattern (or group of patterns) p as the harmonic mean of the precision and recall achieved by the logistic regression classifier, formulated as follows:$${mathrm{FS}}_p = 2 times frac{{({mathrm{precision}}_p times {mathrm{recall}}_p)}}{{({mathrm{precision}}_p + {mathrm{recall}}_p)}}$$
(7)
$${mathrm{Precision}}_p = frac{{{mathrm{TP}}_p}}{{{mathrm{TP}}_p + {mathrm{FP}}_p}}$$
(8)
$${mathrm{Recall}}_p = frac{{{mathrm{TP}}_p}}{{{mathrm{TP}}_p + {mathrm{FN}}_p}}$$
(9)
where, given positive and negative classes, TPp, FPp, and FNp are the number of true positive elements (the classifier correctly predicts the positive class for them), false-positive elements (the classifier incorrectly predicts the positive class), and false negative elements (the classifier incorrectly predicts the negative class). A weighted average of the FSp is then computed to account for class imbalance:$${mathrm{FS}}_{{mathrm{average}}} = frac{1}{H}mathop {sum }limits_{p in P} |p| times {mathrm{FS}}_p$$
(10)
where P is the total number of classes p and H is the total number of elements aggregated across all classes.Software implementationWe coded the exclusion criteria in Matlab and used the “prctile” function for the calculation of the quantiles in Tukey’s fences (last Matlab version tested: R2020b)50. We implemented the customer segmentation analysis and logistic regression classifier in Python (version 3.7.1): the customer segmentation analysis relies on the hierarchical clustering included in the SciPy library51; the logistic regression classifier, along with its k-fold cross-validation and performance evaluation, were implemented using the machine learning library Scikit-learn52; SMOTE oversampling was implemented using the Imbalanced-learn toolbox53. A notebook with the Python code used to generate the results reported in this article is available in a public GitHub repository54. More