Management actions to mitigate the occurrence of pharmaceuticals in river networks in a global change context.

Human consumption of pharmaceuticals leads to high concentrations of pharmaceuticals in wastewater, which is usually not or insufficiently collected and treated before release into freshwater ecosystems. There, pharmaceuticals may pose a threat to aquatic biota. Unfortunately, occurrence data of pharmaceuticals in freshwaters at the global scale is scarce and unevenly distributed, thus preventing the identification of hotspots, the prediction of the impact of Global Change (particularly streamflow and population changes) on their occurrence, and the design of appropriate mitigation actions. Here, we use diclofenac (DCL) as a typical pharmaceutical contaminant, and a global model of DCL chemical fate based on wastewater sanitation, population density and hydrology to estimate current concentrations in the river network, the impact of future changes in runoff and population, and potential mitigation actions in line with the Sustainable Development Goals. Our model is calibrated against measurements available in the literature. We estimate that 2.74 ± 0.63% of global river network length has DCL concentrations exceeding the proposed EU Watch list limit (100 ng L-1). Furthermore, many rivers downstream from highly populated areas show values beyond 1000 ng L-1, particularly those associated to megacities in Asia lacking sufficient wastewater treatment. This situation will worsen with Global Change, as streamflow changes and human population growth will increase the proportion of the river network above 100 ng L-1 up to 3.10 ± 0.72%. Given this background, we assessed feasible source and end-of-pipe mitigation actions, including per capita consumption reduction through eco-directed sustainable prescribing (EDSP), the implementation of the United Nations Sustainable Development Goal (SDG) 6 of halving the proportion of population without access to safely managed sanitation services, and improvement of wastewater treatment plants up to the Swiss standards. Among the considered end-of-pipe mitigation actions, implementation of SDG 6 was the most effective, reducing the proportion of the river network above 100 ng L-1 down to 2.95 ± 0.68%. However, EDSP brought this proportion down to 2.80 ± 0.64%. Overall, our findings indicate that the sole implementation of technological improvements will be insufficient to prevent the expected increase in pharmaceuticals concentration, and that technological solution need to be combined with source mitigation actions.


Introduction
Pharmaceuticals are administered worldwide as prescription medicines, over-the-counter therapeutic drugs, and veterinary drugs. Their consumption and subsequent incomplete assimilation by humans or animals lead to high concentrations in wastewater, which reach freshwater ecosystems either directly as untreated sewerage or open defecation, or after treatment at either in-situ sanitation facilities (e.g. septic tanks) or off-site at wastewater treatment plants (WWTP). Once in the environment, pharmaceuticals might pose a threat to biological communities. Indeed, many pharmaceuticals cause harmful effects on freshwater organisms at environmentally relevant concentrations (Boxall et al., 2012;Brodin et al., 2013). Unfortunately, occurrence data of pharmaceuticals in freshwaters at the global scale is scarce and unevenly distributed (aus der Beek et al., 2016). For example, most occurrence data of the anti-inflammatory DCL (CAS#15307-79-6) are located in Europe and east Asia, but are poor or absent in many other regions (Acuña et al., 2015a,b;Lonappan et al., 2016). Similarly, occurrence data of pharmaceuticals in seawater at the global scale is also scarce and unevenly distributed (Arpin-Pont et al., 2016;Bonnefille et al., 2018). To fill gaps in occurrence data, we need either monitoring surveys or estimating occurrence by means of modeling (Gimeno et al., 2017;Johnson et al., 2013;Ort et al., 2009;Strokal et al., 2019). Because of these knowledge limitations, we cannot anticipate how Global Change (particularly streamflow and population changes) might influence occurrence of pharmaceuticals, and therefore river basin authorities cannot design appropriate mitigation actions.
Here we use DCL as a typical pharmaceutical contaminant and a global model considering human consumption, removal by different sanitation treatment facilities, and natural attenuation in freshwater ecosystems (Font et al., 2019) to estimate DCL concentration in the global river network. We first develop and calibrate the model to estimate DCL concentration under current conditions in order to fill occurrence gaps in the global river network. Second, we use the calibrated model to predict the effects of changing river streamflow and growing human population as a likely Global Change scenario for 2030. Third, we assess the effectiveness of either source or end-of-pipe management actions to mitigate the concentration increases associated with Global Change. Specifically, we assessed the decrease in the per capita consumption of DCL through EDSP, the implementation of the United Nations SDG 6 (UN General Assembly, 2015), and the upgrading of the WWTP to the Swiss current standards for microcontaminants removal (Eggen et al., 2014).

Diclofenac as a model compound
DCL has been selected as a typical pharmaceutical because of the relative abundance of occurrence data and because it is acknowledged as the most commonly used non-steroidal anti-inflammatory drug, with a market share close to that of the next 3 most popular anti-inflammatory drugs combined (ibuprofen, mefenamic acid, and naproxen) (McGettigan and Henry, 2013). Furthermore, anti-inflammatory drugs are among the most prescribed and used pharmaceuticals worldwide because of their efficacy and low abuse potential (Brune and Patrignani, 2015). Finally, DCL is an emerging environmental contaminant included till recently in the EU watch list and on the list of priority substances in the United Kingdom.

Global DCL model
We have used an improved version of the GLOBAL-FATE model (Font et al., 2019), a GIS-based model that simulates the transport and fate of pharmaceuticals for human use in the global river network, considering direct disposal to freshwaters and through wastewater treatment plants. The main improvement over the former version is the inclusion of a third contaminant pathway (in situ treatments) to the river network (see below).

Model structure
Estimating DCL concentration along discretized river networks is based on a local mass balance that relates downstream loads in a river segment or grid cell to loads from upstream segments and local consumption, excretion, and decay in wastewater treatments and the river network, including reservoirs and lakes. The core function of the model is the load calculation: where L j is the DCL load in cell j and N j are the 8 neighbor cells contributing load (L i ) to L j . The contaminant load produced in cell j is a function of population (P j, supplied as a population raster), DCL consumption (M j , supplied at the country level in g inhab -1 y -1 ), and the human excretion rate (supplied as a fraction for the whole Earth). Contaminant load generated by population locally in cell j can experience decay in WWTPs or in situ treatments (including disposal in soils). w wwtp and w rest are the fractions of wastewater treated in WWTP and in situ, respectively, supplied at the country level differentiating between urban and rural areas. wwtp and rest are the pharmaceutical removal rate in WWTP and in other sanitation services, respectively. Finally, the total load in cell j experiences decay in the river network following a first order reaction, controlled by a decay constant (k, hour -1 ) and the residence time of the river water in the cell (RT, hours). The output of the equation is the contaminant load (conveniently transformed from g to ng y -1 ) at each pixel, that is later converted into DCL concentration (ng L -1 ) dividing the load by an estimate of streamflow volume in each cell.
To calculate RT and streamflow at each cell, we routed the runoff generated locally at each cell along the river network. As a runoff raster, we used a published composite global annual runoff (Fekete et al., 2002), which consists in a raster of annual runoff with values in L m −2 year −1 . The runoff was then routed using the flow direction raster of the Dominant River Tracing (DRT) (Wu et al., 2012), a database designed to perform macro scale hydrologic calculations. The streamflow values generated were then used to calculate water velocity solving the Manning equation as in Font et al. (2019), who calculated velocity at each cell deriving slopes from a global digital elevation model, considering a constant Manning coefficient (0.044 s·m −1/3 Schulze et al. (2005)), and estimating the hydraulic radius using the power functions of Leopold and Maddock (1953). RT was estimated from velocity and the length of the river network at each cell. Residence time in lakes and reservoirs were calculated using available global databases on the location, shape, and volume of lakes and reservoirs. These spatially explicit databases are converted into a raster with the same resolution and projection as the other hydrological rasters. Residence time was calculated from volume and streamflow (Font et al., 2019). The general strategy to include lakes in the contaminant simulation is to store all features of a given lake (volume, residence time) in the outlet cell (i.e., the cell routing the streamflow downstream from the lake), making the rest of cells of the lake as mere pipes of water and constituents to that outlet cell, where all contaminant reactions occur (Font et al., 2019).
GLOBAL-FATE spatial resolution for this study was 1/16 degree, corresponding to a cell size of approximately 7 × 7 km at the Equator. For this study, we discarded all arid regions from final calculations, because the runoff database used is highly uncertain in regions showing less than 100 mm yr −1 (Schewe et al., 2014), leading to unreasonable DCL concentrations. Similarly, only DCL concentrations modelled in pixels draining > 150 km 2 (i.e., a pixel receiving two upstream pixels) were kept for final calculations. This is based on the fact that we identified huge DCL concentrations in large urban areas due to unrealistic representation of river reaches and water infrastructure at the working resolution in these areas (sewage infrastructure in large urban areas is not accounted for in our model) (Font et al., 2019). The following is a description of the main data sources used in the implementation of GLOBAL-FATE in this work.

Consumption
We estimate human consumption based on per capita consumption (M) and human population (P). Human population was obtained from the freely downloadable Gridded Population of the World version 4 V. Acuña, et al. Environment International 143 (2020) 105993 (GPWv4) (Center for International Earth Science Information Network -CIESIN -Columbia University, 2016), whereas per capita consumption was estimated from the IMS-Health dataset for the period 2011-2013. Specifically, the IMS-Health dataset includes national human consumption of DCL data for 86 nations (expressed as kilograms of consumed DCL per year). Therefore, we estimated national consumption for the remaining 145 nations. Although IMS-Health data was only available for 38% of the global nations, these included ca. 82% of the global population. National per capita consumption for the 86 nations included in the IMS-Health dataset was estimated as the national consumption divided by the national population. The per capita consumption values of nations not included in the IMS-Health dataset were estimated as equal to the adjacent nation per capita consumption (using Adjacent Fields function of ArcMap, ESRI). DCL is consumed for both human and animal health, and it can be applied on the skin or administered orally in the case of humans. Among these, oral administration it is the main form of administration and accounts for about 70% of the worldwide DCL sales following IMS-Health data (Zhang et al., 2008). The oral administration, either domestic or in hospitals, is the only source of DCL consumption data at the national scale and has been therefore used as input data in the model. However, the IMS-Health database had data gaps, as data for domestic and hospital use was not always reported. In order to correct these gaps, we included 3 correction parameters in the model during calibration. Specifically, the consumption correction parameter 1 (B 1 ) was applied to those countries with data on both domestic and hospital use, so that values above 1 would account for skin application and animal health; the correction parameter 2 (B 2 ) was applied to those countries with only domestic use, and the correction parameter 3 (B 3 ) was applied to those countries with only hospital use data (Table 1).

Removal at sanitation facilities
We used the "Drinking Water Sanitation Hygiene" database from the United Nations World Health Organization (version July 2017) (WHO/UNICEF Joint Monitoring Programme for Water Supply and Sanitation, 2017) to route wastewater from human excretion to the river network through 3 pathways: wastewater transported through a sewer and discharged to the river (i), wastewater transported and treated off-site (w wwtp ) (ii), and the rest (w rest ), which included less than safely managed and treated and disposed in situ (iii). The first was estimated as the difference between the wastewater collected and transported through sewers and the wastewater treated at WWTP. The second was estimated as the sum of wastewater treated at WWTP and "emptied and treated", which accounts for the safely managed wastewater treated off-site. The third was estimated as the sum of "septic tanks", "latrines and other", "unimproved sanitation", "limited (shared) sanitation", and "open defecation". No removal rate was considered for the first route of untreated sewerage, and a percent removal rate was used for the second and third routes ( wwtp and rest ). The initial value and the range of possible removal rates in WWTP (Table 1) to be used during the calibration process were based on a literature review (Tiedeken et al., 2017), whereas the initial value and the range of possible removal attenuation rates in the other sanitation facilities (Table 1) was based on expert knowledge. In fact, the wider range of values reflects the uncertainty associated with this removal rate, which lumps together very different contaminant pathways.

In-stream attenuation
The processes that drive in-stream attenuation (i.e., biological transformations, photolysis, sorption, volatilization) of pharmaceuticals depend on compound characteristics and on several physicochemical and biological parameters including streamflow, temperature, turbidity, dissolved oxygen concentration, biofilm biomass, and pH (Acuña et al., 2015a,b). Consequently, in-stream attenuation rates show high variability consistent with the spatial and temporal heterogeneity of driver variables, which complicates the prediction of attenuation rates across rivers. The paucity of studies reporting in-stream attenuation rates preclude a comprehensive analysis on the variability of this parameter at large scales for a particular compound, such as DCL. Instead, we decided to assign a single in-stream attenuation value for the whole global river network. Although this is a huge simplification at odds with the varying nature of in-stream attenuation rates measured in the field, we decided not to over-parameterize our model considering that the limited information on in-stream concentrations for pharmaceuticals would not allow for a proper identification of different instream attenuation rates depending on potential driver variables (e.g., stream order, streamflow). To identify the best value of the decay constant for first order in-stream decay, we reviewed the literature (Aymerich et al., 2016;Acuña et al., 2015a,b;Johnson et al., 2013;Kunkel and Radke, 2011;Radke et al., 2010) to set an initial value and a range for the calibration (Table 1).

Calibration
The calibration of the model parameters k, wwtp , rest , and the consumption data parameters B B B , , and 1 2 3 (see previous section on consumption correction) was done through the genetic algorithm from the NLOPT library (https://nlopt.readthedocs.io/en/latest/), which delivers a probable distribution of values for each calibrated parameter. The upper and lower bounds for parameters during calibration were decided using previous knowledge as detailed in the preceding sections (Table 1). We used the sum of the absolute differences between the logobserved and log-simulated DCL loads obtained with Eq. (1) as objective function (Acuña et al., 2015a,b). The observed loads were obtained by multiplying observed DCL concentrations and the calculated flows at the corresponding cells. After the calibration process reached the convergence criterion, we assessed the goodness of fit of model predictions using the Nash-Sutcliffe model efficiency coefficient (E), taking the median of the distribution of each calibrated parameter to solve the model ( Table 1). The relationship between observed and predicted values ( Fig. 1) shows an E value of 0.37, which is a reasonably good result. In fact, global models for contaminants always suffer from low to medium performance scores due to the scarce and spatially biased V. Acuña, et al. Environment International 143 (2020) 105993 datasets available for model evaluation , and frequently only show E values > 0.5 after intensive calibration procedures (Harrison et al., 2019). In any case, models of this type aim at reconciling estimates with measurements and at exploring scenarios at large spatial scales, but they cannot be interpreted as simulations of the fate and transport of DCL that mimic reality across scales, particularly at local ones.

Uncertainty analysis
The uncertainty in the GLOBAL-FATE model output has several sources, such as input data and parameters, in addition to model assumptions and simplifications. Considering the complexity of our modeling exercise (11 potential adjustable parameters and 4 global rasters providing key inputs), we discarded the calculation of uncertainties using Monte Carlo techniques. Instead, we used a stochastic approach to estimate the total uncertainty, without separating each individual contribution (Montanari and Brath, 2004). After calibration, the probability distribution of the model error is assessed conditioned to the simulated DCL load, and the uncertainty for a particular prediction is extracted from this probability distribution. In less technical terms, the scatter plot between modeled and observed data (in our case contaminant loads) is used to extract a Gaussian error distribution for any value along the predicted variable. That is, for a particular modeled value, an error distribution is generated using the information contained in the comparison between modeled and observed loads. Considering that our data set is quite extensive (more than 400 data points), representing a wide range of hydrological conditions and geographical locations, the scatter between observed and modeled data is a representation of the error propagated by inefficiencies in all inputs, parameters, and assumptions used in the model. Therefore, the error distributions obtained with this method account for the uncertainty propagated by the whole parametric, input, and model structure. This is a great advantage over Monte Carlo methodologies that tend to restrict uncertainty analyses to adjustable parameters, frequently resulting in unrealistic narrow uncertainty bounds at odds with a visual comparison between modeled and observed values. The drawback is that we cannot identify the main sources of uncertainty in the model outputs. However, considering the tremendous savings in computational time and the wide uncertainty bounds obtained (which place us in a very conservative position respect the conclusions reached with our modeling), we opted for the stochastic approach. We used this method to report all uncertainties related to metrics reported in the paper, i.e. loads, concentrations, and river kilometers with concentration exceeding 30 and 100 ng L -1 .

Limitations of the modeling approach
Considering the spatial resolution of our model (cells are ca. 7 × 7 km in the Equator), and the nature of most of the input information related to contaminant consumption and treatment, we suggest not using results from this implementation of GLOBAL-FATE to draw conclusions at the single cell resolution (Font et al., 2019). Also, GLOBAL-FATE is a steady state model, so it cannot dynamically simulate climatic extreme events or seasonality in hydrology or population. Therefore, all results in this paper refer to average conditions. Our working resolution also implies that our study is blind to impacts on very small rivers, which would be severely affected by sewage discharge contamination. In fact, sewage pollution including DCL in small, low order streams is widespread (Acuña et al., 2015a,b). This would imply that this study is underestimating the actual risk of DCL contamination to surface waters, therefore conclusions in this study can be considered as conservative.
GLOBAL-FATE is among the first contaminant models that fully integrates lakes and reservoirs in the routing of a contaminant along the global river network (Font et al., 2019). The implementation of GLOBAL-FATE in this work also considers the routing of contaminants in lakes and reservoirs, but we could not identify reliable DCL concentration data in our observed data set sampled in lakes, because most samples were collected very close to the shore or to contaminant point sources, thus being a very poor proxy of average DCL concentration in the whole lake, which is what GLOBAL-FATE can simulate. Therefore, while still affecting the routing of contaminants, and in the absence of appropriate observed data to compare with, we do not report results for lakes and reservoirs in this paper. Another major limitation of the study is the focus on a single pharmaceutical compound, with a particular consumption rate per capita and removal efficiencies at sanitation facilities. Thus, other compounds differing in these 2 key aspects would have provided different occurrence patterns, which limits the generalization of our results. However, our focus is not on the simulated spatial patterns, but on the scenario analysis, which should not be dramatically affected by the consumption and attenuation rates. However, we fully support the recent call for multi-pollutant modeling , as the conclusions drawn from a multi-pollutant modeling exercise would be more robust.

Global change scenario
The scenario for our planet in 2030 was defined considering plausible changes in run-off (Schewe et al., 2014), population growth (Jones and O'Neill, 2016), and an estimate of the sanitation services based on the projection of the trends from the years 2000 to 2015 (business-asusual sanitation scenario) (Fig. 2). Note that run-off changes and population growth scenarios were chosen according with the Shared Socioeconomic Pathway 2 (SSP2) (O'Neill et al., 2017), which corresponds to the business-as-usual scenario.

Mitigation scenarios
We simulated the effects of the EDSP as a source mitigation action. Specifically, we implemented a 10% reduction of the per capita consumption of DCL. In fact, DCL can be either sold over-the-counter or after medical prescription, and the balance between these two forms changes over time and across nations (Gimeno et al., 2018). Despite the relevance of the over-the-counter sales in DCL, we applied a reduction of the per capita consumption which has been reported as feasible through EDSP for other pharmaceuticals (Daughton and Ruhoy, 2013). V. Acuña, et al. Environment International 143 (2020) 105993 Another relevant assumption in this mitigation action is that we implemented a standard reduction for all nations instead of implementing a reduction proportional to the current national per capita consumption. Needless to say, modeling results are sensitive to these assumptions, but they constituted the best trade-off between implementing a realistic scenario for pharmaceuticals, and the specific features (i.e. consumption national patterns, and chemical properties) of DCL. Among the SDG of the United Nations, SDG 6 refers to drinking water, sanitation and hygiene. To define our scenario, we only considered targets 6.2 and 6.3, which refer to the access to adequate and equitable sanitation and hygiene for all and ending open defecation (6.2), and to halving the proportion of untreated wastewater and substantially increasing water recycling and safe reuse globally (6.3) (WHO/UNICEF Joint Monitoring Programme for Water Supply and Sanitation, 2017). Given the structure of our model, target 6.2 does not influence the partitioning of wastewater among the 3 wastewater pathways to the river network we model. However, target 6.3 does influence the model reducing by 50% the proportion of untreated wastewater (Fig. 2). This reduction affected all categories other than safely managed, including estimated untreated sewerage. Wastewater was reallocated to either safely treated in-situ or off-site, and the proportion going to each one of these options was different for urban and for rural areas depending on the current proportion of wastewater treated either in situ or off-site. Specifically, in urban and rural areas, all reallocated water from untreated sewerage goes to off-site treatment (i); in rural areas, reallocated water from the third route goes to either in situ (41%) or off-site (59%) (ii); and in urban areas, reallocated water from the third route goes to either in situ (23%) or off-site (77%).
As for possible technological improvements to enhance the removal efficiency in WWTP, we implemented a scenario simulating the Swiss national strategy for upgrading WWTP at the global scale. To identify the WWTP to upgrade, we used the criteria defined at the Swiss strategy: WWTP serving more than 80,000 inhabitants (i), and WWTP serving more than 8000 inhabitants and contributing more than 10% of the dry weather stream flow (ii) (Eggen et al., 2014) (Fig. 2). In those WWTP, the removal efficiency of WWTP was set at 80%, which is a reduction commonly reported when using ozonation followed by sand filtration, an acknowledged technique to eliminate pharmaceuticals from wastewater (Eggen et al., 2014). The cost of implementing this mitigation scenario was estimated multiplying the population affected by the scenario implementation times the cost of treating wastewater with ozonation (1.9 ± 3.1 € inhabitant −1 y -1 ). This value includes construction costs, operational costs (i.e. personnel, maintenance), and variable costs (i.e. electrical consumption for ozone generation and sand filtration, as well as the cost of pure oxygen for the ozone production). Specifically, cost per inhabitant was estimated as the product of wastewater generation per inhabitant per year times ozonation treatment cost per volume of water, which is inversely proportional to the population served by the WWTP (Gimeno et al., 2018). This assumes an ozone dosage of 0.7 g O 3 per gram of dissolved organic carbon, and a retention time in the ozonation tank of 25 min.

Results and discussion
Our results indicate that 2119 Mg of DCL are consumed by the World population every year, a number higher than previous estimates (1450 Mg y -1 for 2012, (Acuña et al., 2015a,b); and 940 Mg y -1 for 2005 (Zhang et al., 2008)). Our estimate differed from previously reported ones because we allowed the model to correct for data inhomogeneities in the data-bases used (see methods). 264 Mg y -1 are excreted and introduced into the global water cycle, of which 30.4 Mg y -1 are removed in WWTP, 11.9 Mg y -1 are removed by in-situ treatment facilities or by natural soils when disposed in situ (see Section 2.2.4.), 88.8 Mg y -1 are attenuated in the river network, and the remaining 133 Mg y -1 are delivered to the oceans. The estimated global WWTP removal efficiency was 40.4% (Table 1), a value slightly lower than the mean but within the range reported for conventional treatment facilities (63 ± 24%) (Tiedeken et al., 2017) and doubling a previous estimate based on  V. Acuña, et al. Environment International 143 (2020) 105993 Europe (Johnson et al., 2013). National DCL removal by sanitation services ranged from 0 to 40% and was closely related with the per capita nominal GDP (Fig. 3a). Thus, relatively low removals by sanitation services occur in nations under development, and removal values above 30% are only observed in developed nations (Supplementary  Table S1). In contrast with the relatively high removal efficiency at WWTP, the estimated removal efficiency of in-situ treatment facilities or by natural soils when disposed in situ was only 8.4% (Table 1), indicating that transport and treatment of wastewater in WWTP involves a major change in the attenuation of pharmaceuticals occurrence in freshwaters. The estimated natural attenuation in rivers (k = 0.00032 d -1 ; Supplementary Table S1) was lower than those reported in studies at the river segment scale (Aymerich et al., 2016), although higher than others that assumed DCL to behave conservatively (Johnson et al., 2013;Kunkel and Radke, 2011;Radke et al., 2010). Therefore, although the attenuation rate in freshwaters is small, the V. Acuña, et al. Environment International 143 (2020) 105993 amount of DCL attenuated at the global scale is by far larger than that removed by sanitation services, and it is therefore advised to always consider attenuation by freshwater ecosystems when modeling DCL chemical fate. The median DCL concentration in the global river network is 0.27 ng L -1 , and ranges from 0 to 233 ng L -1 (99% distribution bounds) (results available online in raster format at https://github.com/icra/ DCL ( Fig. 4a and Supplementary Table S2). Concentration is particularly high in extensive areas of the Indian subcontinent, East Asia, Central Europe, and Tropical Africa. In fact, many rivers downstream from highly populated areas show values beyond 1000 ng L -1 ; particularly those associated with megacities in Asia without sufficient wastewater treatment. This is also reflected at the river basin scale, as Asian basins like the Yellow River or the Ganges have the highest mean concentrations, although European river basins also show high values (Table 2). Overall, our results indicate that around 2.7 ± 0.6% of the global river network (i.e., 172,114 km) has DCL concentrations exceeding the proposed EU Watch list limit (100 ng L -1 ) (Council of the European Communities, 2013). Likewise, 3.9 ± 0.9% of the EU river network (i.e., 12,293 km) has DCL concentrations exceeding 100 ng L -1 , similar to previous reports for the European Union, stating values between > 1 and 10% (Johnson et al., 2013). Moreover, 7.7 ± 0.9% of the global river network exceeds 30 ng L -1 , which is the concentration identified as a threshold for measurable toxic effects in freshwater organisms (Acuña et al., 2015a,b). Last but not least, we estimated that 22.9% of the global coastline exports some DCL. Specifically, the mass of DCL delivered to the oceans differs considerably among and within countries, with a median of 1.03 g per coastline km and year, and ranging from 0 to 2,474 g km −1 y -1 (99% distribution bounds) (Fig. 4b).
Differences also exist between river basins, with values usually below 1 Mg y -1 (Table 2), but with basins such as the Ganges exporting around 5% of the total global export. These exports are responsible of the DCL occurrence in oceans, which range from a few ng L -1 to 100 ng L -1 (Bonnefille et al., 2018).
The modeled Global Change scenario for 2030, which considered the joined effects of changes in run-off (Schewe et al., 2014), changes in population (Jones and O'Neill, 2016), and an estimate of the sanitation services for 2030 based on the projection of trends from 2000 to 2015 (business-as-usual sanitation scenario, see details in supplementary information), predicted an increase of the proportion of the river network exceeding 100 ng L -1 up to 3.10 ± 0.72%. Large differences will occur between continents, as the increases in concentration mainly occur in sub-Saharan Africa, the Mediterranean, the Caribbean, and some southeastern countries (Fig. 4c), while in other regions DCL concentration will even decrease, due to higher dilution capacity with higher streamflow values. These differences are also reflected at the river basin scale, with some basins experiencing increases above 30% and others decreases below 15% (Table 2). Global Change will also considerably increase the exports to oceans in most river basins (Table 2).
Decreasing the per capita consumption through EDSP brought the proportion of river network exceeding 100 ng L -1 down to 2.80 ± 0.64% (Supplementary Table S2). Regarding the end-of-pipe mitigation actions, the implementation of SDG 6 brought the proportion of river network exceeding 100 ng L -1 down to 2.95 ± 0.68%. The implementation of tertiary treatment according to Swiss standards in WWTP affected 13.7% of the global population, and brought the proportion of river network exceeding 100 ng L -1 down to 2.98 ± 0.62%. Finally, the combined effect of all considered mitigation actions resulted in a proportion of the global river network exceeding 100 ng L -1 of 2.54 ± 0.58%; which implies that the effects of Global Change in the occurrence of DCL were only effectively mitigated when implementing all the mitigation actions together. When analyzing the results at the scale of SDG regions (WHO/UNICEF Joint Monitoring Programme for Water Supply and Sanitation, 2017), EDSP is the most effective mitigation action except in Eastern Asia and South-eastern Asia, and in Australia and New Zealand, (Fig. 5). In regards to end-ofpipe mitigation actions, the implementation of SDG6 was considerably more effective than the Swiss standards in the regions under development (SDG regions Sub-Saharan Africa, Western Asia and Northern Africa, and Central Asia and Southern Asia) (Fig. 5), where the implementation of advanced tertiary treatments for microcontaminants removal is often limited due to lack of existing WWTP.
The implementation costs differ considerably between mitigation actions. The estimated cost of implementing the SDG 6 scenario is 65.2 · 10 9 € y -1 (Hutton and Varughese, 2016), whereas implementing the Swiss standards globally would cost 5.5 · 10 9 € y -1 . Therefore, the implementation of the Swiss standards was more cost-effective than the implementation of the SDG 6 in mitigating the effects of Global Change, although it should be stressed that this comparison is only for a single organic microcontaminant, and that a different pattern might emerge  V. Acuña, et al. Environment International 143 (2020) 105993 when considering all potential contaminants. Interestingly, the implementation cost of the Swiss standards in WWTP is directly proportional to nominal GDP (Fig. 3b). Thus, the effort for each nation would be proportional to their income, facilitating this significant improvement in sanitation services.

Conclusions
Overall, the sole implementation of technological improvements will not suffice to prevent the expected increases of the concentration of pharmaceuticals such as DCL in many regions, and affordable source mitigation scenarios such as the simulated EDSP, environmental friendly drug design and process development, and the take-back and management of unused drugs (Daughton and Ruhoy, 2013) need to be considered. Another reason for combining source and end-of-pipe mitigation actions is that the projected impacts of Global Change are largest in developing regions (Fig. 4c), where wastewater treatment is still limited and where it might be difficult to deal with the costs associated to fully implementing the SDG6 (Guiteras et al., 2015;Kaiser, 2015;Mara and Evans, 2017;Narain, 2010). Similar calls for coordinated efforts have been already done (Wen et al., 2017), and here we note that the combination of actions seems to be the only reliable strategy to counteract the Global Change-led impact of pharmaceuticals on global freshwaters.