Environmental risk assessment of fluctuating diazinon concentrations in an urban and agricultural catchment using toxicokinetic-toxicodynamic modeling.

Temporally resolved environmental risk assessment of fluctuating concentrations of micropollutants is presented. We separated the prediction of toxicity over time from the extrapolation from one to many species and from acute to sublethal effects. A toxicokinetic-toxicodynamic (TKTD) model predicted toxicity caused by fluctuating concentrations of diazinon, measured by time-resolved sampling over 108 days from three locations in a stream network, representing urban, agricultural and mixed land use. We calculated extrapolation factors to quantify variation in toxicity among species and effect types based on available toxicity data, while correcting for different test durations with the TKTD model. Sampling from the distribution of extrapolation factors and prediction of time-resolved toxicity with the TKTD model facilitated subsequent calculation of the risk of undesired toxic events. Approximately one-fifth of aquatic organisms were at risk and fluctuating concentrations were more toxic than their averages. Contribution of urban and agricultural sources of diazinon to the overall risk varied. Thus using fixed concentrations as water quality criteria appears overly simplistic because it ignores the temporal dimension of toxicity. However, the improved prediction of toxicity for fluctuating concentrations may be small compared to uncertainty due to limited diversity of toxicity data to base the extrapolation factors on.

Introduction

Background

Concentrations of micropollutants(1) such as pesticides and biocides in freshwater streams may strongly fluctuate over time, as demonstrated by an increasing number of chemical monitoring studies with high temporal resolution[2−8] or predicted by pollutant fate models.(9) Assessing the risk of adverse effects of such fluctuating concentrations to aquatic organisms is challenging,[2,7,10−13] because traditional risk assessment methods are based on ecotoxicological tests, employing constant exposure concentrations and fixed durations, that do not explicitly consider the temporal aspect of toxicity. As a result, the current risk assessment procedures lack methods to assess the potential toxicity of fluctuating and repeated pulsed exposures.

Traditionally derived water quality criteria for short-term or long-term exposure do not suffice because it is unclear how to evaluate very short peaks above the short-term quality criterion or concentrations between the short-term and the long-term criterion.[2,7,14] In addition, it is unclear how to deal with subsequent pollutant peaks, which raises the question whether exposed organisms have had enough time to recover in between exposures,[7,14,15] and whether low concentrations between peaks contribute to the risk of adverse effects.[7,14]

Risk assessment of fluctuating concentrations requires assessment of toxicity over time (prediction of toxic effects beyond test conditions), and such a prediction requires a mechanistic model of toxicity.[10,16] Toxicokinetic-toxicodynamic (TKTD) models (14,17) are suitable for this purpose as they can simulate time-dependent phenomena related to prediction of toxic effects such as carry-over toxicity and cumulative effects,[15,16,18−23] delayed and post exposure effects (19,24−27) and organism recovery.[18,28−31]

Objectives

The objective of this study was to overcome the limitations of time-invariable toxicity data by implementing simulated time-variable toxicity into risk assessment of micropollutants. We demonstrate that prediction of toxic effects over time is feasible with a mechanistic model being used for two purposes: (i) the simulation of effects over time and (ii) the calculation of extrapolation factors for a set of ecotoxicological data. We used the insecticide diazinon as an example and assessed the contribution of urban and agricultural sources of diazinon toward the ecotoxicological risk based on time-resolved exposure data, measured at different locations within a catchment.

Outline of the Time-Resolved Risk Assessment Approach

The goal of the time-resolved risk assessment is to predict the risk posed to aquatic organisms as a function of cumulative exposure over time. We illustrate the approach using measured time series of diazinon exposure concentrations from three locations in a small stream network, which represent urban, agricultural and mixed land use, as inputs for the environmental risk assessment.

Risk assessment of micropollutants encompasses three dimensions: (i) prediction of toxicity under long-term and fluctuating exposure, (ii) extrapolation from lethal to sublethal end points and (iii) extrapolation from one (or a few) tested species to many (all) species in the water body of concern. Here, we separated the three dimensions and addressed the prediction of toxicity over time separately from the extrapolation to many species and sublethal effects. Prediction of toxic effects over time was achieved with simulations by a mechanistic TKTD model.

Fixed assessment factors are traditionally used to address the above-mentioned extrapolation steps in risk assessment of chemicals.[32−36] Here, we replace these assessment factors with extrapolation factors that quantify the variation in the extrapolation of toxicity from one to many species and from lethal to sublethal effects based on available, relevant toxicity data. In the calculation of these extrapolation factors, we corrected effect concentrations for differences in test durations by using the mechanistic TKTD model. Finally, we fitted a distribution to the extrapolation factors, sampled a large number of extrapolation factors from that distribution and predicted the time course of toxic effects with the TKTD model for each of these extrapolation factors. We then calculated the risk of undesired toxic events from the multitude of these simulations (Figure 1).

Figure 1

Calculation of the effect and risk curves (middle and upper panels in Figure 2). Toxic effects from time variable exposure patterns are predicted with a mechanistic effect model (TKTD model). Extrapolation between species and from lethal to sublethal effects is achieved by multiplying the exposure concentration (C(t)) with extrapolation factors. The extrapolation factors account for the variation in sensitivity of different species and effect types, based on the data in Table 1 (see Scheme 1 for explanation of the calculation steps).

Figure 2

Measured concentrations of diazinon (lower panels) in different parts of the catchment, effect curves simulated for different percentiles (solid lines from left to right: 95th, 90th, 87.5th, 85th, and 80th percentile) of the extrapolation factor distribution (middle panels) and the risk (top panels). Risk is the fraction of probabilistic simulations that show toxic effects, i.e. the fraction of affected combinations of species and effect type (see Figure 1 and Scheme 1 for calculation steps). In the lower panels, a peak originating from the agricultural part of the catchment is indicated by (a) and elevated concentrations originating from the urban part by (b). Note the different scale of the y-axis in the lower panel (station AGR). Day 0 corresponds to the March 10, 2007. The effect curve of the example AA-EQS value is plotted in the middle panels (dashed line).

Table 1

Ecotoxicological Studies Relevant in the Context of the Protection Goals of the Swiss Water Protection Law(39)a

species	effect type	effect measure	effect size (% affected)	reference for effect size	test duration (days)	effect concentration (nmol/L)	ref	corresponding LCx for same effect size and duration in G. pulexnmol/L	extrapolation factor (EF)	log EFc
Brachionus calyciflorus	reproduction	NOEC	10.5	(42)	2	26 286	RIVM database (41)	96.5	0.0037	–2.340
Brachionus calyciflorus	reproduction	EC50	50	EC50	2	36 144	RIVM database (41)	262.3	0.0073
Brachionus calyciflorus	reproduction	MATC	25	estimatedb	2	32 858	RIVM database (41)	150.1	0.0046
Ceriodaphnia dubia	reproduction	NOEC	25	(42)	7	0.72	RIVM database (41)	24.2	33.48	1.421
Ceriodaphnia dubia	reproduction	EC43	43	(61)	7	1.31	(61)	26.8	20.39
Ceriodaphnia dubia	reproduction	EC28	28	(61)	7	0.95	(61)	25.1	26.34
Daphnia magna	filtration rate	EC50	50	EC50	1	1.54	RIVM database (41)	1279.4	828.5	2.918
Daphnia magna	reproduction	NOEC	10	(42)	21	0.66	RIVM database (41)	8.9	13.48	1.130
Daphnia pulex	reproductive rate, population extinction	EC99	99	(62)	2	6.57	(62)	1289.4	196.2	2.293
Daphnia pulex	population size	EC09	9	(63)	10	2.04	(63)	13.7	6.725	0.828
Pimephales promelas	growth	NOEC	10.5	(42)	7	282.9	RIVM database (41)	19.3	0.0681	–1.228
Pimephales promelas	growth	NOEC	10.5	(42)	32	164.3	RIVM database (41)	7.8	0.0474
Pimephales promelas	growth	NOEC	10.5	(42)	32	131.4	RIVM database (41)	7.8	0.0592
Pimephales promelas	reproduction	NOEC	14	(43)	168	10.51	RIVM database (41)	7.0	0.6657	–0.177
Pogonichthys macrolepidotus	spinal deformities	EC35	35	(64)	4	4107	(64)	58.2	0.0142	–1.849
Oncorynchus mykiss	swimming	EC40	40	(65)	1	1643	(65)	999.2	0.6082	–0.216
Salvelinus fontinalis	reproduction	NOEC	14	(43)	11	1.81	RIVM database (41)	13.6	7.548	0.878

We calculate the concentration that causes the same percentage mortality in G. pulex (LCx) for the same duration and percentage effect as in each of the ecotoxicological studies using the TKTD model. Extrapolation factors are the ratio of the LCx and the effect concentration in the ecotoxicological study. Multiple values of the same combination of species and effect type are replaced with their median.

As the MATC is the geometric mean between NOEC and LOEC a larger effect size than the NOEC can be expected, thus we use the largest effect size from all NOEC end points here, i.e. 25%.

A normal distribution with mean = 0.3325 and SD = 1.667 was fitted to the log EFs in the last column.

Scheme 1

Step by Step Explanation of the Risk Assessment Calculations

Step (1) is the calibration of the toxicokinetic-toxicodynamic (TKTD) model.(19) Step (2) generates the data for the probabilistic risk assessment, which is carried out in steps (3) to (6). The TKTD model is used in steps (3) and (5).

Materials and Methods

Diazinon Monitoring Data

Diazinon is a hydrophobic insecticide (log Kow 3.81(37)), which frequently occurs in pulses or fluctuating concentrations in streams[3,5,6,8] and its metabolite diazoxon inhibits acetylcholinesterase. Next to its use in agriculture, diazinon is also used as insecticide in urban areas (e.g., against lice on roses, fish moths in wet rooms etc.). The diazinon data originate from the study described in Wittmer et al. (2010).(8) The study area (Figure S1 in Supporting Information (SI)) is a subcatchment embedded in the catchment of Lake Greifensee located on the Swiss Plateau. The studied catchment covers 25 km2 and has a mixed urban and agricultural land use. We compared the temporal dynamics of diazinon concentrations in a subcatchment with only agricultural land use (AGR in Figure S1 of the SI) to one with predominant urban land use (URB in Figure S1 of the SI, including effluent from the wastewater treatment plant (WWTP)). Settlements located in the subcatchment AGR are connected to the sewer system in URB. There is no urban stormwater discharged into the river in AGR. We compared the dynamics of the two subcatchments to the dynamics at the outlet of the entire catchment (TOT in Figure S1 of the SI), which comprises 470 ha of arable land and 12 000 inhabitants (10 000 in the subcatchment URB).

The three stations were monitored from 10th of March to 26th of June 2007. During several rain events throughout the monitoring period, samples were taken at a high resolution (15 min to hourly composite samples). Weekly grab samples were taken during dry weather conditions. Samples were stored at minus 20 °C in the dark to prevent degradation. Prior to analysis an isotope-labeled internal standard was added to 50 mL of filtered sample. Analytes were enriched with online solid phase extraction, separated by liquid chromatography, and detected by tandem mass spectrometry (SPE-LC-MS/MS). For details on the analytical method, see Singer et al. (2010)(38) and for details on sampling procedure and sites see Wittmer et al. (2010).(8)

Ecotoxicological Data

We demonstrate the method on the example of a catchment in Switzerland and Swiss protection goals. Swiss protection goals for surface water aim to protect fisheries, the health of animals and biological processes needed by plants and animals to fulfill their physiological needs, specifically metabolism, reproduction, and olfactory orientation.[39,40] Thus we considered sublethal toxic effects as relevant and corresponding to the protection goals and collected such data for diazinon and aquatic organisms as input for this case study. Lethal data were not used in calculations, because assessment based on sublethal data is also protective for effects on survival. Risk assessments for different protection goals may use different criteria to select input data.

Ecotoxicological data on sublethal effects of diazinon in aquatic organisms were collected from two sources. The RIVM database contains data from the U.S. EPA and The Netherlands National Institute of Public Health and the Environment (RIVM)(41) until the year 2000 (D. de Zwart Pers. comm.) and is quality controlled (see de Zwart(41) for procedure). From this database we selected all entries for sublethal effects on freshwater organisms (i.e., excluding mortality, immobility, abundance, and undefined effect types), which yielded twelve records. From these records, we deleted data on algae (1 record) because they do not carry the target enzyme, the acetylcholinesterase, and are not as sensitive as fish and invertebrates toward diazinon. Effect sizes for NOEC values were taken from Crane and Newman(42) and Suter et al.(43) Additional data were collected from peer-reviewed literature by searching the SCOPUS database. The query for “diazinon” “toxicity” “*water” AND (year after 2000) in title, abstract, or keywords yielded 102 abstracts (15 February 2011). From these, we selected studies with information relevant to the Swiss water protection goals yielding six additional records. Here, effect sizes were taken from the original references. All studies are listed in Table 1.

Toxicokinetic-Toxicodynamic Model

The TKTD model was previously parametrized to simulate survival of the stream dwelling arthropod Gammarus pulex.(19) The freshwater arthropod G. pulex is frequently used in ecotoxicological studies,(44) is very sensitive to diazinon (96 h-LC50 is 4.15 μg/L(45)) and exhibits similar sensitivity to environmental pollutants as Daphnia magna.(46)G. pulex is abundant in headwater streams of our case study catchment. Furthermore, its longevity makes it a suitable test organism to parametrize a TKTD model, which captures the relationship between time and toxicity for diazinon, including carry-over toxicity.(19)G. pulex requires approximately 28 days to recover from toxic stress caused by diazinon.(19) If G. pulex is pre-exposed to diazinon, then a subsequent exposure can cause stronger toxic effects than what would be expected from a one–pulse exposure alone.(19) Such carry–over toxicity was caused by slow organism recovery as shown by the TKTD model.(19) The TKTD model for diazinon and G. pulex simulates the time-course of toxicant uptake, biotransformation and elimination as a first step (toxicokinetics, TK) and the development of damage within the organisms and subsequent increased mortality as a second step (toxicodynamics, TD).(19) This TKTD model corresponds to the special case of stochastic death in the General Unified Threshold Model of Survival.(17) The model was parametrized using measured internal concentrations of diazinon and its biotransformation products, in particular diazoxon, within G. pulex as well as long-term survival experiments with pulsed exposure patterns.(19) This mechanistic model describes the processes leading to toxicity on a temporal scale and so captures the time-toxicity relationship for survival of G. pulex and diazinon (step 1 in Scheme 1). Thus, the model can be used to predict survival of G. pulex for other, untested diazinon exposure patterns. Here we also use the TKTD model as a proxy for the time-toxicity relationships of other combinations of species and effect type, because the TKTD model for survival of G. pulex is the best approximation for these unknown time-toxicity relationships that we have, even if that requires some bold assumptions.

We used the TKTD model for two purposes. First, we calculate extrapolation factors for extrapolation from G. pulex to other species and effect types. Second, we simulate the time-course of risk of adverse toxic effects for aquatic organisms. Details of both calculations follow below.

Assumptions about the Relationship between Exposure Time and Toxicity

Assumptions about the relationship between exposure time and toxicity are rarely made explicit in current risk assessment schemes. The default method for dealing with time-variable exposure concentrations in risk assessment is to use time-weighted average concentrations.[10,47,48] This method is based on Haber’s law and assumes that the product of exposure time and concentration determines toxicity, i.e., the same time integral of concentration yields the same toxicity,[10,49] and that carry-over toxicity does not occur. Consequently, the average concentration of a fluctuating exposure pattern is assumed to result in the same toxicity as the fluctuating exposure. Deviations from Haber’s law can be caused by toxicokinetics (carry-over bioaccumulation), toxicodynamics (carry-over toxicity) or both. Relationships between time and toxicity can be modeled with TKTD models,(17) such as the one we base the proposed risk assessment on.(19) We assume that the relationship between time and toxicity in all the species that we address with our risk assessment is the same as that in G. pulex. This does not mean that we assume the same sensitivity, but we assume the same organism recovery time.(18) Assuming that all assessed species resemble G. pulex in their time-toxicity relationship is the best currently possible approximation, because to date G. pulex is the only species for which a time toxicity relationship of diazinon was quantified with a TKTD model. We further assume that an extrapolation factor that we multiply exposure concentrations with, can account for sensitivity variation of different species and effect types. This is necessary because we lack the knowledge to incorporate possible nonlinearity in that aspect of sensitivity variation. Note that applying species sensitivity distributions or assessment factors to averaged exposure concentrations, as in standard procedures,[35,47,48] also assumes the same time-toxicity relationship for all organisms and it ignores the possibility of carry-over toxicity.

Extrapolation Factors

We define extrapolation factors (EFs) as ratio between the lethal concentration for x% of G. pulex at duration t (LCx(t)) and the effect concentration for x% effect in the ecotoxicological studies with different species and effect types at duration t (ECx(t)) (Table 1, step 3 in Scheme 1). Multiplying the effect concentration ECx(t)i from a study in Table 1 with the extrapolation factor EFi, results in the corresponding lethal concentration LCx(t)i for G. pulex for the same duration t and same effect size (x%).

Derivation of EFs for the ecotoxicological studies in Table 1 (step 2 in Scheme 1) requires first calculating the LCx(t)i of G. pulex that corresponds to the same duration t and x percentage effect of the sublethal effects ECx(t)i on other species than G. pulex. This calculation was carried out using the TKTD model: survival of G. pulex under a constant exposure concentration LCx(t)i was simulated and the concentration adjusted such that the effect level on mortality of G. pulex corresponded to the effect level (x %) of the sublethal effects ECx(t)i on the other species after the duration t. EFi was calculated as ratio LCx(t)i/ECx(t)i.

All EFs were log transformed to achieve normally distributed data (passed Kolmogorov–Smirnov test, D’Agostino and Pearson omnibus normality test and Shapiro-Wilk normality test). Multiple values of EFs for the same combination of species and effect type were replaced with their median (last column in Table 1). The log-transformed EFs in our case study are characterized by a normal distribution with mean 0.3325 and standard deviation 1.667 (n = 11). This is step 4 in Scheme 1.

We also calculated EFs for mortality caused by diazinon in different aquatic species in the RIVM database and compared their distribution to sublethal EFs (see SI). Here, only the distribution of sublethal EFs was used in further analyses (see rationale in section Ecotoxicological data).

With the assumptions that the time-effect relationship for G. pulex and diazinon resembles that of the other combinations of species and effect type and that the EF captures sensitivity variation, one can do the following calculations: First multiply the concentration time series of diazinon in any exposure pattern with the EF for a given combination of species and effect type. Then use that concentration time series as input in the TKTD model calibrated for G. pulex survival and simulate the time-course of effects for the combination of species and effect type of that EF.

Effect Simulations

We repeatedly simulated survival of G. pulex in response to measured concentration time series in a small stream multiplied with the EF. Without multiplying exposure with the EF the simulation would predict survival of G. pulex. With the EF the simulated effect curve is the prediction for a sublethal effect on another species. If the simulation were done with one of the EFs from Table 1, then the simulated effect curve would be the time course of effect for that specific combination of biological species and sublethal effect type. The application of the extrapolation factor changed the interpretation of the model output from simulated time-course of survival of G. pulex into simulated toxic effects for the respective combination of aquatic organisms and effect type.

We randomly sampled 1000 log EFs from the normal distribution of log-transformed sublethal EFs, back-transformed them to EFs (denoted EF), multiplied the concentrations in the exposure time series with these EFs and simulated the time-course of effect with the TKTD model for each of these instances (see Scheme 1, step 5). The same 1000 EFs were used for each of the three concentration time-series corresponding to the total catchment and its urban and agricultural subcatchments, respectively.

The survival probability S(t), which is the output of the original TKTD model,(19) was then converted to effect as follows:where EFn [−] is the extrapolation factor in simulation run n (out of 1000 runs), C(t) [nmol/L] is the measured concentration time series of diazinon(8) and Effect [−] is the predicted effect size. Effect size has the same interpretation as in the ecotoxicological study that corresponds to the EF, e.g., at the EF of 0.6082 a simulated effect of 0.1 could be interpreted as swimming being affected in 10% of Oncorhynchus mykiss (see Table 1).

Calculation of Risk

Risk can be defined as the probability of an undesired event.(50) Here we assume that the ecotoxicological data in Table 1 represents variability in sublethal toxicity of diazinon to fish and aquatic invertebrate species that the Swiss protection goal aims to protect. Thus we define an undesired event as one where Effect(t) exceeds 50%, i.e., where 50% of the organisms of a species suffer from an undesired toxic effect (one could also chose a different effect size). As we simulated a sample of 1000 effect curves based on the distribution of EFs, we can calculate risk as the fraction of those simulations where Effect(t) exceeds 50% (Scheme 1, step 6). As each of the 1000 EFs drawn from the distribution represents a single, possible combination of species and type of sublethal effect, this is the fraction of affected combinations of species and effect type, i.e., risk (top panels in Figure 2).

Results and Discussion

The toxicity data from literature and thereof derived EFs are listed in Table 1. The distributions of EFs for lethal and sublethal effects are plotted in Figure S1 of the SI. Both distributions span about 5 orders of magnitude. Even though EFs for sublethal effects are larger by approximately 1 order of magnitude, a considerable overlap in the two distributions is obvious. This comparison indicates that differences among species contribute more to variation than differences between lethal and sublethal effects.

Time-Course of Effect and Risk to Aquatic Organisms

The measured diazinon concentrations fluctuated strongly (Figure 2, lower panels). At the outlet of the catchment (station TOT, Figure S2 of the SI), elevated background concentrations between 0.5 and 1 nmol/L were observed throughout most of the study with pronounced peaks occurring at irregular intervals. Discharge of station URB (Figure S2 of the SI) consists mostly of urban stormwater runoff and effluent of the WWTP. There was no high peak at the beginning of the time series, as in the agricultural part of the catchment, but background concentrations were higher than at the outlet of the entire catchment (TOT) because there was less dilution of the WWTP discharge by discharge from agricultural areas.

Using measured time-courses of diazinon concentrations as input for modeling resulted in a risk of adverse effects to 17% (TOT), 20% (URB) and 24% (AGR) of aquatic organisms at the three monitoring stations (Figure 2, top panels). Interestingly, the temporal onset and contribution of urban and agricultural sources to the risk varied between sites.

Selected percentiles of the effect curves predicted using the TKTD model are plotted in Figure 2, middle panels. Effect curves show the increasing effects over time for the chosen percentiles of combinations of species and effect type in relation to the concentration curve plotted below. From left to right (total, urban, agricultural part of the catchment) the exposure profiles affect a larger fraction of combinations of species and effect type and more severe effects are seen earlier. In the agricultural exposure profile the first high peak is the dominant toxic event, whereas in the other two profiles toxicity results from combination of lower long-term exposure and short peaks.

Low long-term and low-level pollution from urban sources of diazinon, which was observed at stations TOT and URB, caused risk for the 5% most sensitive combinations of species and sublethal effect type (Figure 2, middle panels, 95th percentile curve). The high peak around day 8, coming from the agricultural part of the catchment (labeled (a) in Figure 2, lower panels), raised the risk in station TOT to 11%. The lower exposure levels between days 25 and 45 increased the risk by another 2% and then peaks around day 60 raised the risk to the final total of 17% in station TOT. After day 50 risk at station URB increased to higher percentiles than at station TOT, i.e., toward the end of the monitoring period, we predict that more combinations of species and effect type are affected at station URB (20%) than station TOT (17%).

A different situation was observed at the monitoring station AGR, where discharge consists only of agricultural and some road runoff. Hardly any constant background concentrations were observed. Instead short concentration peaks occurred (Figure 2, lower panels). Being about 1 order of magnitude higher than at the other two stations, the short pulse (1 to 2d) around day 8, caused most likely by agricultural activities alone, led to 24% of risk.

Risk Assessment of Time-Varying Exposure

The concept of using fixed concentrations as water quality criteria is overly simplistic because it ignores the temporal dimension of toxicity that is evident from the examples above. Risk assessment based on TKTD models demonstrates that the risk of adverse effects can be caused by short, high peak exposures or by lower, long-term exposures or by a combination of both. Due to the mechanistic nature of the underlying TKTD model, which has been shown to properly predict carry-over toxicity,(19) we are confident that the risk assessment procedure presented here captures the temporal aspects of toxicity. Alternative methods (reviewed elsewhere[2,7,10]), such as time-weighted averages, are less well qualified for risk assessment of fluctuating concentrations, because they are not able to predict carry-over toxicity.(16) In this case study, the average concentration underestimated toxicity by up to a factor of 4 (see below). Thus a dynamic, mechanistic effect model complements classical risk assessment based on short- and long-term quality criteria.

Quantification of spatial and temporal exposure can be carried out together with toxicokinetic–toxicodynamic modeling of effects, which places the temporal extrapolation of toxicity on a mechanistic basis. The approach presented here separates the temporal extrapolation step from the extrapolation to other biological species and other end points (e.g., sublethal effects). Consequently, the separated extrapolation steps can be tested, re-evaluated, improved, and validated independently of each other.

The use of the EF for probabilistic risk assessment is in analogy to the species sensitivity distribution approach.(51) As with species sensitivity distributions, the available relevant toxicity data determines the outcome of the risk assessment. The toxicity data (Table 1) as input for the risk assessment model, can be easily exchanged by a different set of toxicity data, perhaps to reflect a different protection goal and risk assessment context. The combination of TKTD models and probabilistic calculations allows quantification of the time course of effect curves for different combinations of species and effect type. Further, it allows calculation of the time course of risk in response to fluctuating exposure concentrations. The proposed time-resolved risk assessment method is generic and can be adopted for risk assessment of other pesticides(2) and biocides,(52) or the assessment of pollution in the context of the Water Framework Directive.(53)

Comparison with Effects Based on Average Concentrations

We quantified the difference in toxicity of the fluctuating exposure patterns compared with the toxicity of their average concentrations. For each of the three different exposure patterns and for each of the three corresponding average exposure concentrations, we calculated the factor that the exposure time series must be multiplied with to result in 50% mortality of G. pulex at the end.

The concentration time series from the total catchment had to be multiplied with 75 to result in 50% mortality at the end (107.4 d), whereas the corresponding average concentration of 0.0760 nmol/L required multiplication by a factor of 96 for the same toxic effect. For the fluctuating exposure time series from the urban part of the catchment the factor was 46 and for the corresponding average concentration of 0.1187 nmol/L the factor was 62. For the fluctuating exposure time series from the agricultural part of the catchment the factor was 26 and for the corresponding average concentration of 0.0726 nmol/L the factor was 101. Thus, in all three exposure patterns, the fluctuating concentrations were more toxic than their averages (1.28-, 1.35-, and 3.88-fold, respectively). Here, the TKTD model was used to simulate toxicity for the full length of the exposure profile. In practice, the time-weighted average concentrations are calculated for a shorter period of time(2) and how large the resulting error of these methods is needs to be systematically investigated.

Comparison with Environmental Quality Standards

For comparison with environmental quality standards derived with the commonly applied method in the European Union,(54) we took the Annual Average Environmental Quality Standard (AA-EQS) of 0.015 μg/L (0.049 nmol/L) from a report(55) that compared different methods to derive environmental quality standards for Swiss surface waters. Assuming 365 days as a test duration and an effect size of 10% for the AA-EQS, which is generally derived from NOEC data, we calculated the corresponding EF as 141. This EF corresponds to the 88th percentile of the EF distribution, thus it protects 88% of combinations of species and effect type based on the data in Table 1. The corresponding effect curves are plotted in the middle panels in Figure 2 (AA-EQS, dotted line). The average concentrations in stations TOT, URB and AGR were 0.076 nmol/L, 0.119 nmol/L and 0.073 nmol/L, respectively. Thus comparison with the AA-EQS of 0.049 nmol/L already indicates a risk. The added value of the probabilistic assessment is the quantification of total risk from any exposure pattern and the possibility to investigate effects on specified percentiles of species effect combinations.

Uncertainty and Limitations

Ecotoxicological data for this case study comprises only three taxonomic groups: rotifers (3 studies), crustaceans (7 studies), and fish (7 studies). We do not know how well this limited diversity of taxa in our data captures variation in sensitivity between different species. Similarly diversity of effect types in our data set is limited and biased toward reproduction in rotifers and crustaceans and growth in fish. As with any modeling study, quality of the results depends on the input data. In this case, the limited diversity of species and effect types in the ecotoxicity database poses the greatest source of uncertainty. The desired diversity in ecotoxicity data was also discussed in the context of species sensitivity distributions(51) and these discussions suggest, that also for studies like the one presented here, a higher diversity of ecotoxicological data is needed.

Variation in sensitivity of different species and effect types spans several orders of magnitude, whereas the difference between the toxicity of the fluctuating exposure profiles and that of their corresponding average concentrations was less than 4-fold. Thus we improved prediction of toxicity for fluctuating concentrations with our method, but these improvements may be small compared to the uncertainty due to limited diversity of toxicity data. Consequently, the potential of this new method can best be realized in combination with larger sets of ecotoxicity data. For small sets of ecotoxicity data large uncertainties remain.

Relevant Sources and Priorities for Mitigation Measures

Urban and agricultural sources of diazinon contributed to the risk of adverse effects on aquatic organisms in the studied catchment. Their respective contribution varied and cannot be quantified a priori or as any fixed ratio. Rather, both sources contributed with variable proportions and the relationship by which long-term pollution and short concentration peaks induce adverse effects is nonlinear. From a practical point of view, the analysis demonstrated that, depending on the location in the river network, either urban or agricultural sources of diazinon or their interplay caused a risk of adverse effects. Without the TKTD model, just relying on the observed concentration data, it was difficult to judge the ecotoxicological relevance of the different exposure patterns and sources. Through the TKTD modeling, we learn that an integrated approach to risk mitigation that targets both urban and agricultural sources is most likely to efficiently reduce risk of adverse effects to aquatic organisms. Hence it is important to understand the dynamics of discharge from different sources.(8) For diazinon, restrictions or mitigation measures on both urban and agricultural uses, are required.

Spatial Aspects of Risk Assessment in Aquatic Systems

Availability of measured concentrations time-series at different locations in the catchment raises the question how to aggregate risks calculated at various locations in any given river network. Our set of three locations illustrates that the same event, for example the high peak originating from the agricultural area around day 8, may cause strong risk of adverse effects in one location (24% in AGR) but much less in another (11% in TOT). Integration of TKTD models with spatially explicit models for populations, such as individual or agent based models,[56,57] may be one possible solution if the assessment is carried out at the population level. If however, the protection goal is aimed at the organism level, then the question of how to aggregate risk of adverse effects from different locations remains.

Research Needs

TKTD models can handle fluctuating concentrations of micropollutants. However, the approach is currently limited by the availability of model parameters. These parameters are unique for each combination of chemical and biological species and more of these relationships of time and toxicity need to be quantified. The current reliance on one such relationship for all other assessed species carries uncertainty, which we need to quantify and reduce. If TKTD models for more species, compounds, and effect types were available, then these could ultimately replace the extrapolation factors used here. Raw data from standard acute toxicity tests can already be used for parametrization of TKTD models(17) and selection of species and sublethal end points for TKTD model development(58) could be prioritized based on ecological considerations (e.g., trophic position, ecosystem services) to maximize relevance.

We assessed the risk from one compound alone, although it is known that multiple chemical stressors occur simultaneously or after each other, also in the catchment studied here.(8) Risk assessment of mixtures in time is important as different compounds can interact and contribute to mixture toxicity even if they occur days(22) or weeks apart.(15) TKTD models for mixtures of chemicals are available,[15,59] but their relationship to standard models for toxicity of mixtures with simultaneous, constant exposure is not well established.(60) TKTD models for temporal mixtures can explain phenomena such as the sequence effect,(15) which is caused by carry-over toxicity, however the large number of model parameters needed to assess the temporal mixtures of chemicals on a catchment scale are not available yet.