Spatial and temporal analyses of perfluorooctanoic acid in drinking water for external exposure assessment in the Ruhr metropolitan area, Germany

Perfluorooctanoic acid (PFOA) and related chemicals among the per- and polyfluoroalkyl substances are widely distributed in the environment. Adverse health effects may occur even at low exposure levels. A large-scale contamination of drinking water resources, especially the rivers Möhne and Ruhr, was detected in North Rhine-Westphalia, Germany, in summer 2006. As a result, concentration data are available from the water supply stations along these rivers and partly from the water network of areas supplied by them. Measurements started after the contamination’s discovery. In addition, there are sparse data from stations in other regions. Further information on the supply structure (river system, station-to-area relations) and expert statements on contamination risks are available. Within the first state-wide environmental-epidemiological study on the general population, these data are temporally and spatially modelled to assign estimated exposure values to the resident population. A generalized linear model with an inverse link offers consistent temporal approaches to model each station’s PFOA data along the river Ruhr and copes with a steeply decreasing temporal data pattern at mainly affected locations. The river’s segments between the main junctions are the most important factor to explain the spatial structure, besides local effects. Deductions from supply stations to areas and, therefore, to the residents’ risk are possible via estimated supply proportions. The resulting potential correlation structure of the supply areas is dominated by the common water supply from the Ruhr. Other areas are often isolated and, therefore, need to be modelled separately. The contamination is homogeneous within most of the areas.


Toxicological background
Per-and polyfluoroalkyl substances (PFASs) name a group of man-made persistent organic chemicals produced since the 1940s which are used in a wide variety of products used by consumers and industry (OECD 2004;Buck et al. 2011). Their biochemical persistence and stability have led to pollution of the environment worldwide and an internal exposure of the general human population (cf., e.g., Buck et al. 2011;OECD 2004). Perfluorooctanoic acid (PFOA) and perfluorooctane sulphonic acid (PFOS) are readily absorbed after ingestion or inhalation, not metabolized and Electronic supplementary material The online version of this article (https://doi.org/10.1007/s00477-020-01932-8) contains supplementary material, which is available to authorized users. slowly excreted via urine and feces. Mean human half-lives are estimated to be 3. 5-3.8 and 4.8-5.4 years, for PFOA and PFOS, respectively (Olsen et al. 2007). The human population's most important source is food, with contaminated drinking water predominating if present. Therefore, an increased water contamination may be used as a surrogate marker of the human internal exposure, where this is considerably above background level: Consumption of contaminated drinking water is associated with increased PFOA-concentrations in human plasma . Tap water contributes to more than 75% of the estimated daily intake in comparison to other exposure pathways (diet, air, consumer articles) in highly affected regions (Vestergren and Cousins 2009).
The toxicity of PFASs is currently being discussed among toxicologists, epidemiologists and regulatory agencies. Developmental and liver toxicity (Lau et al. 2004;Lau 2012), elevated blood lipids (Frisbee et al. 2010;Steenland et al. 2009) and immunotoxicity (NTP 2016) seem to be major health effects associated with exposure to PFOA and PFOS. Other relevant impairments to health, e.g. the incidence of type 2 diabetes in humans (Sun et al. 2018), have been reported. Meta-analyses of animal experimental and human epidemiological data concluded that developmental exposure to PFOA adversely affects human health based on sufficient evidence of decreased fetal growth in both human and non-human mammalian species Lam et al. 2014). However, the mechanisms of toxicity have not yet been sufficiently clarified. Detailed toxicological assessments have been published lately (EPA 2016; ATSDR 2018; EFSA CONTAM 2018; UBA HBM Commission 2018; NJDEP 2019).
The European Food Safety Authority (EFSA) lowered the tolerable daily intake from 1500 (150) (EFSA 2008) to 0.9 (1.9) ng per kg body weight for PFOA (PFOS) and concluded that a considerable proportion of the population currently is exceeding these levels (EFSA CONTAM 2018). Only recently, EFSA also evaluated the exposure to the sum of four different PFASs and recommended an even lower tolerable weekly intake of 4.4 ng per kg body weight for the sum of PFOA, PFOS, PFHxS and PFNA (EFSA -CONTAM 2020). The German Drinking Water Commission lowered their action guidance value for pregnant women, breastfeeding mothers and infants under 24 months of age to 50 ng PFOA per litre of drinking water (Trinkwasserkommission 2020). Among the effects evaluated by the EFSA and the drinking water commission are lipid metabolism, immunotoxicity and developmental effects.

Contamination event and drinking water data base
A large-scale contamination of drinking water with PFOA and other PFASs has been discovered in the summer of 2006 in North Rhine-Westphalia (NRW), Germany (Skutlarek et al. , 2006). The principal cause has been the use of PFASs-polluted conditioner (fertilizer) on more than 1300 farmlands with subsequent contamination of surface waters. Further PFOA-contaminations of both surface and groundwater by sewage plants and other sources like fire extinguishing foam have occurred in different parts of NRW (LANUV 2011), though with minor effects on drinking water. PFASs are generally not removed during riverbank filtration and conventional water processing; some water supply stations have subsequently installed activated charcoal filters in order to reduce the PFOA-concentrations in drinking water below the guidance value of the German Drinking Water Commission (Trinkwasserkommission 2016).
Since 2006, the NRW state environmental agency has implemented an extensive PFASs monitoring programme for relevant environmental media, including soil, water and drinking water (LANUV 2011). The resulting drinking water database has been complemented with data collected by the authors from water suppliers. These data constitute the basis for our analyses.
A particularly severe contamination has affected the river Ruhr and its tributary Möhne (Skutlarek et al. 2006;Wilhelm et al. 2008)). Both rivers are important resources of the regional drinking water: More than 20 important water supply stations are installed along them, supplying a large region of NRW, especially the Ruhr metropolitan area with about a quarter of NRW's 18 million inhabitants.
Human biomonitoring studies have been conducted (Hölzer et al. , 2009) with a focus on the town of Arnsberg, supplied by a Möhne-dependent water supply station near the river's mouth, where PFOA concentrations above 500 ng/l have been measured. (The lifelong tolerable  guide  value  is  currently  set  to 100 ng/l, (Trinkwasserkommission 2016). Furthermore, adverse health effects are suspected to occur even in low ranges of PFOA exposure.)

Aims
This article describes the first phase of our PerSpat (Perfluoroalkyl Spatial) project, which has been designed to assess the exposure of the residential population to PFOA, in order to analyse possible associations with relevant health indicators like birth registry data by means of temporal and spatial modelling in a second phase.
This environmental-epidemiological study on the general population, living in a large region of variable and, in parts, very high PFOA concentrations in drinking water, is the first of this kind on state level and, therefore, importantly contributes to the epidemiological research on developmental toxicity of PFASs.
This article contains the first comprehensive description of the collected and processed complex secondary data and accompanying information on PFOA concentrations measured in drinking water following the NRW contamination incident. In this analysis, we focus on PFOA among the PFASs due to higher data variability, measurements considerably above critical values, and its long half-life.
We evaluate statistical models to find adequate representations of these measurements and to predict the PFOA concentration at particular water supply stations and points of time. Thereby, we encounter a new and unusual data situation in terms of the combination of very instationary temporal processes (in terms of both mean and variance), two spatial levels of interest, and the connection via rivers.
After fitting Generalized Linear Models (GLMs) for temporal regression at the stations, we explore the spatial structure of the results along the well-observed river Ruhr.
A second spatial data level-drinking water supply areas-is distinguished. We numerically estimate the proportions of drinking water from different stations supplying the same area, using meta information on the NRW water supply structure, in order to assess the regional average PFOA concentrations.

Related modelling approaches
We start our empirical data analysis rather open-minded with regard to model classes. Ultimately, we aim at a model in continuous time and discrete space.
As for approaches in similar scenarios and profound model developments, see, e.g., White et al. (2017) for spatial smoothing models for areas being but partially observed: Data are primarily given on area level with a Markov property assumed, whereas we have to focus on the temporal modelling first and also encounter complicated spatial relationships. We consider approaches like these in future stages of analysis.
Bayesian hierarchical discrete-space approaches, along with intrinsic CAR prior choice considerations, can be found in Keefe et al. (2019). Again, the neighbourhood structure is more specific than found in our situation. Grollemund et al. (2019) present a regression model using functional predictors. These could be representations of the temporal structures in a future spatio-temporal analysis in our situation; before, we have to explore the temporal data patters.
As a frequentistic alternative, Tang et al. (2019) develop a semi-parametric copula model to cope with the spatiotemporal dependence and to predict values for new points of time and sites. This is rather restrictive, being Markovian in time, whereas we model continuous, irregularly measured temporal data. Wang and Sun (2019) propose a spatial regression model with local polynomials modelling the coefficients as a function of space, but without an explicit temporal component.
The set-up of consistent spatial correlation matrices along more complex river networks than ours is described as a particular challenge (cf., e.g., Peterson et al. 2007). With more data than we have from groundwater-stations not directly adjacent to a river, the latter could be regarded as an inhomogeneous line source of contamination (cf., e.g., Ayub et al. 2019). Another option would be the incorporation of more detailed river network features, such as water amounts, stream velocity and smaller tributaries (cf., Ver Hoef et al. 2006), if corresponding data were available.
Further methods and application scenarios for risk assessment, when measurement data are available at more steps along water supply systems than in our situation, can be found in Roozbahani et al. (2013). Causal remote effect estimation, when both influence and outcome quantities are measured along a river and over time, unlike in our situation, is performed by Saul et al. (2019) via the parametric g-formula.
Article overview Below, the data are described in more detail in Sect. 2, with an overview on their amount and features (Sect. 2.1), an outline of the NRW river network and its usage for drinking water (Sect. 2.2) as well as the typical temporal structures to be found in the various locations (Sect. 2.3). In Sect. 3, the data are explored with regard to both time and space: Regression models on time are presented in Sect. 3.1, with some findings concerning the spatial structure along the river Ruhr. Section 3.2 deals with the question whether and how to deduce from water supply stations to the contamination of their supplied areas. Section 4 concludes with a discussion, including some preliminary spatial modelling approaches, and perspectives.

Data description
The PFOA measurements have been conducted by the local drinking water suppliers; the data have been collected by the NRW State Agency for Nature, Environment and Consumer Protection (Landesamt fu¨r Natur, Umwelt und Verbraucherschutz; LANUV) and by the Department of Environmental Medicine of the Ruhr-University Bochum.
All data preprocessing steps, visualisations and analyses have been conducted using the R environment (R Core Team 2020) along with its basic packages. Processing of geographical data and drawing of maps have been conducted using the R packages rgdal (Bivand et al. 2019) and tmap (Tennekes 2018).
Some of the maps have been produced using geographical data from the German Federal Agency for Cartography and Geodesy (Ó GeoBasis-DE / BKG, 2019; license: www.gov data.de/dl-de/by-2-0; d a t a : wms_vg250-ew, wfs_dlm1000. Geographical data on the NRW drinking water supply (Figs. 1 and 2) have been provided by the LANUV.

Data amount and features
The data consist of measurements from completely purified, ready-to-drink water in both water supply stations (3349 observations) and networks (536) since summer 2006. Observations are available until 2016 and partially ongoing. Furthermore, there are some samples drawn from raw water (418) and during water treatment (87) at the stations. Many samples, especially the latter, are non-detects. See Fig. 2 for a cross-sectional overview.
NRW is geographically divided in some 450 water supply areas, with any data available from about 200. For the stations, many of them being very small, any data are available from about 250 of 700. There is a complex, neither injective nor surjective assignment between areas and stations, with some minor changes in areal division, station running and assignment over the years. The number of nonequidistant measurements per time period varies between locations and periods according to the contamination risk.
More detailed information on the amount of data per station and their ranges and temporal patterns can be found in Table 1. According to this, the data availability is satisfying for all affected stations, i.e., along the river Ruhr (a total of 2484 observations at all of the 27 stations), and all of them feature at least rather high values; stations not depending on the river Ruhr are much less observed, typically with low values, although there are exceptions. Table 2 details the data availability with respect to water supply areas: Where some of the stations supplying a  particular area are Ruhr dependent, there are typically enough data to assess the area's contamination risk by deduction from the stations. An important issue with regard to modelling are outstanding values at the beginning of the observed period of time: The most affected supply station of Möhnebogen at the Ruhr's tributary Möhne near Arnsberg features a PFOA maximum of 640 ng/l. On the other hand, values below the limit of quantification (LoQ, mostly 10 ng/l) are frequent in less affected regions and later periods of time. A certain discretisation is caused by a data accuracy of 10 ng/l frequently given.
Apart from the measurement data, there is further information from water supply operators on the most likely non-existence of PFOA contamination at certain stations or areas over longer periods of time. For some places, where there are but a few values, or even just one, the measurement usually confirms the water supplier's assumption of no relevant PFASs contamination.
With data being on two spatial levels, stations and areas, there is also information on their assignment and the supplied and demanded water amounts.
In summary, the PFOA concentrations of all water supply stations along the river Ruhr are well observed, beside water network data from many water supply areas depending on them. For other regions, there are fewer measurements per station-usually non-detects-or none. The relationship of stations and areas is explored in Sect. 3.2.1.

Spatial connectivities via the river network
The main rivers of NRW, with water supply stations along them, are shown in Fig. 1. (A graph of their directions of flow and their junctions is shown in the Online Resource, Fig. 1: a simplified spatial structure, comprising of river segments. As explored in Sects. 3.1 and 3.2.2, a certain homogeneity may be presumed for river-depending stations within the same segment.) For the river Ruhr, there are five water supply stations upstream beyond the Möhne mouth, the station of Möhnebogen on the river Möhne, and twenty-one stations downstream from both rivers' junction. With the main PFOA source drained off by the Möhne, the former five stations are less affected and may be neglected in some analyses, yielding a one-dimensional spatial structure, if needed for simplified models.
Apart from the main rivers, there are many stations taking water from smaller tributaries. A weak dependence between these and the respective main river segments may be presumed.

Temporal data patterns for individual water supply stations
We focus on water supply stations along the river Ruhr, with numerous measurements covering a longer period of time. For data availability and patterns, cf., Table 1. According to this, there are two main types of temporal data patterns, a rather diffuse one (Fig. 3) and a striking decrease (Fig. 5).
The first typical temporal structure is found at stations less affected by the contamination incident, being upstream beyond the Möhne mouth (like Mengesohl, Fig. 3 left), partly supplied by groundwater (like Volmarstein, centre) or far downstream at the river Ruhr (like Styrum West, right). Measurements observed here feature no or hardly detectable trends and comparatively small variability. A certain decline in background pollution may be reflected, too.
A second type is a more or less smooth decline over the years, though some short-term trends may be distinguished. The decrease is considerably stronger for the first couple of years (e.g., Essen-Ü berruhr, Fig. 5).
This decrease at the beginning is the strongest at the outstanding station of Mo¨hnebogen, where some water filtering measures have taken place after discovery of the contamination. Filters eliminate all PFASs after installation, become gradually less efficient and have to be reactivated after a while (Fig. 4). Thus, some outliers, change points, and even seasonality may be included in analyses of Table 2 Number of water supply areas by contamination risk and data availability per station: Considering the water supply stations with a relevant supply for a particular area (estimated share of more than 5%, Sect. 3.2.1), we distinguish whether all or some of these stations are potentially affected by the Möhne-Ruhr contamination incident, and whether any data are available from all or some of these stations (right the respective numbers of areas with only one supplying station) this station-and of some others downstream which have installed filters too, if filtering data were available. Apart from these differing observed types of temporal structures, there are periods of missing data in some locations: The density of measurements in the course of time and space varies widely. Thus, extrapolation should be performed with caution, especially when using rather simple regression or interpolation methods.
In summary, there are two typical trends in the PFOA values of water supply stations at the river Ruhr-beside the special case of Möhnebogen station. In particular, a characteristic form of 'decreasing decrease' can be distinguished. Apart from the Ruhr-depending stations, there are but a few with sufficient data to identify a trend (one example are the stationary measurements around the limit of quantification in Haltern, see Fig. 9).

Data exploration
The first part of this Section deals with the temporal modelling of the PFOA concentration data per water supply station along the river Ruhr: The suitability of various regression models is studied and the prediction results of the most appropriate one are discussed with regard to spatial structures along the river. The water supply areas are considered in Sect. 3.2, especially whether and how statements on their PFOA contamination can be deduced from their respective supplying stations.
Let X it denote the PFOA concentration at water supply station i ¼ 1; . . .; m (in our situation, m % 700) for time t ¼ 1; . . .; T. Reasonable units of time are months or days. Let j ¼ 1; . . .; n (n % 450) denote the water supply area. The data are continuous with respect to time with an irregular measurement process: There are periods of time with a high 'data density' (up to two measurements within a few days) and others without any data for several month. Therefore, we use regression models with the time as a covariate instead of autoregressive models. We fit models individually to the data of a water supply station i, but are mainly interested in obtaining a consistent approach along the river Ruhr.
An appropriate model has to cope with the varying shapes of decrease (cf., Figs. 3 and 5), which may also imply some heteroscedasticity. We consider several Generalized Linear Models (GLMs , Table 3) with various distribution families for X it and link functions as well as data transformations.
For the latter, a pre-transformation like lnðX it Þ, ffiffiffiffiffi ffi X it p or 1 ffiffiffiffi X it p and usage of a simple linear regression leads to very large residuals for periods of high measurements. For the pre-transformation 1 X it , the instability of the estimation is even worse. Additionally, results are difficult to interpret, when re-transformed after the model fit.
Instead of that, while sticking with the simple linear regression, we experimentally transform the time covariate: Among the usual fractional polynomials, the model For GLMs without pre-transformation of data, at first we identify all combinations of continuous distribution families and link functions such that the model fit converges for all stations (shown in Table 3). Among these, a Normal distribution with an inverse link turns out to yield the best fit in terms of the residual sum of squares for most of the stations, especially those with the typical 'decreasing decrease'. But also the Gamma distribution with its natural inverse link is often close to it in that sense. Table 3 shows the goodness of fit of the various models. Fig. 5 illustrates their suitability with regard to the estimated mean function as well as their capability to reflect the data's variance. While the variance VarðlÞ of the mean response estimator is always a function of time and is further affected in case of non-identity link functions, a proper choice of the data variance function turns out to be more relevant in terms of scale: VarðX it Þ is constant for all l in case of the Normal family, leading to a certain probability of negative predictions, but also to a closer prediction interval in the period of high values. For the Gamma family, it holds VarðX it Þ ¼ / 2 l 2 ðtÞ, so we are able to interpret and estimate the measurements' standard error as a constant multiple of the mean, in accordance with some assumptions on measurement technique and natural fluctuation (LANUV 2011); on the other hand, the prediction interval is very large in the period of high values (perhaps corresponding with inferior measurement technique in the early period) and rather small later on (notwithstanding that an increased uncertainty is assumed for measurements close the limit of quantification). However, all models sufficiently cover the data variance. The inverse link is very suitable for modelling the 'decreasing decrease' type of temporal data patterns; for less definite trends, especially along the upper Ruhr, others would do (cf. Table 3). In some cases, a slight heteroscedasticity remains. Outlying data at the measurements' beginning have a strong leverage that may lead to a certain instability, especially when 'predicting the past'. With the inverse link being superior compared to the log link, the decrease turns out to be steeper than exponential for highly contaminated places.
Any such models are applicable for situations without change points or outliers as described in Sect. 2.3. Where carbon filtering interventions are known, a piece-wise regression is applied.
With regard to easy interpretability of the parameters, the goodness of the predictions, and the perspective model extension (say, by additional spatial terms), we focus on the Normal GLM with an inverse link. Figure 6 shows a cross-sectionlðt Ã Þ of the mean response estimation results along the river Ruhr. (In the Online Resource, further points of time can be found Fig. 2, a longitudinal overview in Fig. 3 at the same place.)

Discussion of regression results
The results show a clear distinction between the river's segments: The abrupt increase of PFOA concentrations below the junction of Ruhr and Möhne is reflected as well as a decrease below the mouth of the Ruhr's tributary Lenne. The estimations for most of the stations along the lower Ruhr tend to feature smaller uncertainty intervals.
Within the segments, no definite trend is recognizable; local effects of the individual stations seem to overlay a trend along the river. Potential causes may be the varying parts of groundwater in relation to surface water used by the water supply stations as well as their filtering measures. Another aspect are further minor PFOA discharges, especially from sewage plants affecting the river Ruhr (cf. LANUV 2011). In fact, there seems to be a slightly increasing trend along parts of the lower Ruhr.
To conclude, the GLM with Normal distribution and an inverse link function is a convenient way of functional representation of temporal trends in PFOA data of the water supply stations along the river Ruhr. This also opens a possibility for integrating spatial trends in the model. The latter should respect the river segments.

Station-to-area supply proportions
With data and models on two spatial levels-water supply stations and water supply areas-we are strongly interested in the proportions, how much of an area's water stems from the respective stations, in order to mediate between the two levels. However, there is no direct information on this question-apart from rough estimations by a few local drinking water suppliers. Instead of this, we estimate these supply proportions using data on the stations' water supply amounts and the areas' water demands together with the qualitative information which stations supply which areas.
Let a ¼ a 1 ; . . .; a m ð Þ T be the supply and b ¼ b 1 ; . . .; b n ð Þ T the demands. The data (say, cubic metres per year) are somewhat vague and, therefore, standardized to ; respectively, in order to obtain matching totals. Formally, we are in search of a weight matrix such that area j gets the w ji 's part of its water from station i, which means that the water amounts w ji b j of the areas sum up to the stations total supply with the area-wise restrictions Such a deconvolution problem (cf., e.g., Cutler 1978) is generally not uniquely solvable. In our situation, the problem is simplified by several facts: • Only a small number of stations and areas are connected at all. And these relations are known at least qualitatively, albeit not in numbers. Therefore, it is known for far the most of the cases, that a station i Ã does not supply an area j Ã and we have many further restrictions w j Ã i Ã ¼ 0 when calculating W. • If an area is supplied by a single station, the respective weight w j Ã i Ã ¼ 1 is known prior to solving the system of (1) and (2), even if this relation is not bijective.
• Due to the many zeros in W, the equation system of (1) and (2) is actually split in, often very small, partial problems solving the non-zero blocks W ðlÞ 2 0; 1 ½ n l Âm l of W concerning m l stations and n l areas. Moreover, a large number of supply relations are unique, i.e., m l ¼ n l ¼ 1 (cf., Table 4).
Transforming W to a diagonal structure of non-zero blocks -with the help of the R packages igraph (Csardi and Nepusz 2006) and Matrix (Bates and Mächler 2019)we obtain L ¼ 322 linear equation systems like (3). Of these, 27 are non-trivial, i.e., m l ; n l ! 2. Each comprises of m l þ n l independent equations to solve a certain number of unknown quantities w ji 2 ð0; 1Þ.
In 6 cases, the system is determined. In 20 systems, there are more equations than unknowns and we can find a numerical solution for the least squares problem a ðlÞ À W T ðlÞ b ðlÞ 2 2 À! min W ðlÞ !; W ðlÞ 1 m l ¼ 1 n l using the R package NlcOptim (Chen and Yin 2019) to perform optimisation with constraints; the obtained estimated proportions are exact enough given the vague data; in particular, the by far largest system (m l ¼ 72, n l ¼ 82), containing many of the Ruhr dependent stations, is solvable. Only one system is under-determined: it concerns a small region (m l ¼ 3, n l ¼ 6) not related to the Ruhr and is not considered further. The obtained proportions of water supply from stations to areas are important to estimate the stations' relevance for the NRW water supply in general. They are also used to assess the areas' PFOA burden and thereby the residents' risk by weighted averaging. This is particularly relevant if the PFOA concentration is inhomogeneous within an area, especially one supplied from both the Ruhr and other rivers or groundwater (Sect. 3.2.2).
The supply matrix W also illustrates to what extent the spatial units (water supply areas) are contiguous in terms of shared water: Fig. 7 shows one large region of potentially correlated areas (with the river Ruhr in their south) and several smaller clusters. On the other hand, there is a large number of single isolated areas. Therefore, the obtained supply matrix is a key to understand the unusual discrete spatial structure of NRWs drinking water supply. It is the basis for models on supply area level, including analyses of network data.

Data homogeneity of water supply areas
In many cases, PFOA values for water supply areas have to be deduced from the relevant water supply stations' data in order to assess the residents' exposure. We consider whether the involved stations' temporal trends are similar enough ('homogeneous') for joint modelling and whether these data and the respective network data, if any, are sufficiently homogeneous in the course of time.
Homogeneity or inhomogeneity are presumed on the basis of the exploratory results from above, particularly on the two distinguishable main types of temporal data patterns at the stations (Sect. 2.3), the stations' locations with respect to the affected river Ruhr (cf., Sect. 2.2) and the structures found in predictions for Ruhr-dependent stations (Sect. 3.1.2). Table 1, the occurrence of the typical data patterns is strongly correlated with the stations' dependence from the river Ruhr, where high values with a steep decrease in the course of time can be found. Moreover, these patterns are most similar for stations within the same segment of the river, regarding both data and predictions, where prediction intervals often overlap or are near to each other, along the course of time, as shown in Fig. 6 and in the Online Resource, Fig. 2. In other locations, there are typically low values and few measurements. Therefore, we presume an area as homogeneous, if it is essentially supplied by stations, which are in the same situation with regard to their location and their temporal data pattern. Examples for the homogeneity of data, especially of such with a decreasing pattern, are shown in Fig. 8.

As illustrated in
On the other hand, there are a few areas with some supplying stations depending on the contaminated river Ruhr and some being non-affected (cf., the middle part of Table 2). An example of an inhomogeneously supplied area is shown in Fig. 9: Three stations supply the COE_3 area, two of them are dependent on the contaminated river  ( [ 5%) shares of water from the same station): A coloured area is correlated with at least one other area of the same colour and not correlated with any area of other colours. White coloured areas are not correlated with any other area Ruhr, the third (with a water share of about 40%) is nonaffected.
When reliably double-checking the homogeneity, we are limited to cases with sufficient amounts of data from the respective stations or from the area's water network, representing the same periods of time. Larger amounts of station data are mostly given along the river Ruhr (cf., Tables 1 and 2). Network data are rare, especially for the most interesting period of time. (E.g., there are none for the interesting case of the inhomogeneous COE_3 area.) Where available, the range of the network data seems not to substantially differ from the respective concurrent station samples (cf., Fig. 8).
For Ruhr-dependent areas, we combine data from the respective stations, to search for clusters in the space of time and values (using the R package EMCluster, Chen and Maitra 2015). The same is done with network data, where available. It is not possible to distinguish the data and (re-)discover the stations in this way.
Another approach to check whether concurrent station data are distinguishable is a regression model similar to those in Sect. 3.1 for the combined measurement data of an area's supplying stations, with additional dummy variables representing the stations. For areas completely dependent on (a section of) the river Ruhr, the resulting station effects are not significant or at most very weak. For inhomogeneous areas, we can find such effects.
Both approaches let us conclude for the vast majority of water supply areas (cf., the upper and lower part of Table 2), that there is no evidence against a hypothesis of homogeneity, where adequate data are available (i.e., some of the about 30 stations with many data according to Table 1 are involved). This avoids uncertainty by averaging PFOA concentrations. Perspectively, we are able to fit joint models to the respective water supply stations, respecting those without any data, to obtain estimates for a whole area's exposure. Where waters from the river Ruhr and other sources are involved (about 9 areas according to Table 2), this situation is more difficult.

Discussion and conclusions
Given the potential health risks of PFASs, especially PFOA, and the contamination incident affecting the rivers Möhne and Ruhr, we explore the measurement data from water supply stations and within the drinking water network. The data comprise of quite complete series along the river and only few measurements each at other water supply stations. For the Ruhr, we observe a steeply decreasing temporal data pattern at mainly affected stations and a more diffuse structure upstream; later after the contamination, the decrease is weaker. Elsewhere, there are usually non-detects.
We are interested in a consistent temporal modelling approach for comparison of results along the river Ruhr and later building of joint models. A Gaussian Generalized Linear Model with an inverse link function turns out to have the best fit to the data among regression approaches, hinting to a more than exponential contamination decrease after the incident. If there are crucial events such like change points or extreme outliers, a piece-wise regression is applied.
When exploring the spatial structure of the predictions from this model, we find the important role of the river segments (between the main junctions): The PFOA concentration considerably varies between the Ruhr segments, but less and without clear trends within, where station-related effects and minor additional PFOA discharges may prevail.
We are ultimately interested in the residents' PFOA risk assessment and, therefore, have to deduce from station level data and predictions to the water supply areas. We use an approximate supply proportion matrix obtained from rough information of the water amounts. Furthermore, from our exploratory results, we can presume homogeneity of data from stations along the river Ruhr, supplying the same areas, and find no evidence against it. Therefore, the uncertainty when aggregating data in such areas is low. For those supplied by both the river Ruhr and unaffected stations, it is higher, but aggregation is possible using weights from the proportion matrix.
The potential state-wide spatial correlation structure is explored using the supply proportions. We find one large region of potentially correlated areas near the river Ruhr, but also many other supply areas which are isolated or grouped to small clusters. This information is important to set up a covariance matrix for usage in spatial models, but there are still relevant additional dependencies, induced by river or groundwater connections, which are not directly included.
Useful predictions are not possible for largely unobserved regions or periods of time. Extrapolations are challenges for all methods, especially for the early period prior to the first measurements in 2006. Yet, the duration of the Brilon contamination may be estimated to several years from the fertilizer's quantity, the drainage, and the  Fig. 9 Measurement series of PFOA of the three water supply stations supplying the COE_3 water supply area residents' internal exposure (Hölzer et al. 2009;Skutlarek et al. 2006). A lack of data in vast regions also affects the applicability of comprehensive state-wide smoothing models. As a preliminary spatio-temporal approach, we have used a Bayesian conjugate-prior Gamma-Gamma model focussing on the mean surface to predict PFOA values for any water supply station, water supply area and point of time; its uncertainty varies widely according to the presence of data. So, such models are restricted to temporal analyses and to stations along the well-observed river Ruhr. Moreover, water supply areas which are isolated within the supply relationship structure and additional spatial dependencies would have to be respected.
Kriging approaches-using the R packages Com-pRandFld (Padoan and Bevilacqua 2015) and geoR (Ribeiro et al. 2020)-along the river turn out to be inappropriate, too, since they do not respect the river segments, which dominate the spatial data structure. Furthermore, even for the well observed river Ruhr, the data are too sparse to analyse proper variograms, and sometimes the assumption of isotropy seems violated. Ordinary and universal Kriging lead to sufficient local estimates when crossvalidated along the river Ruhr, but perform worse near the Möhnebogen water supply station. It is questionable whether the continuous Kriging solution is an efficient way to merely estimate concentrations at discrete, definite stations, as in our situation.
Topological (Top-) Kriging (using the R package rtop, Skøien et al. 2014;cf., e.g., Laaha et al. 2014) is designed for prediction on the level of the rivers' catchment areas. This would be useful to model measurements along the NRW river network, thereby summarizing several stations each. Nonetheless, with the stations' values connected to water supply areal units, it does not seem productive to introduce a third spatial level. Furthermore, we have no geographical data on the actual catchment areas.
In general and state-wide, the problems of a varying data density, up to entirely non-observed regions, and of a regionally different strength of dependence remain hindering for all Kriging approaches. For the river Ruhr, the rather complicated models, especially Top-Kriging, seem exaggerated for an 'linear' spatial river structure. Therefore, we drop these approaches because of the non-suitable spatial data structure. For problems of the performance of spatial methods, especially (Top-)Kriging, compared to regression approaches for river-related catchment area data, see, e.g., Brunner et al. (2018).
Our station-wise regression results lead to the perspective to combine the temporal regression models for the Ruhr data and add a spatial structure. Spatial effects of the river segments and of the individual water supply stations-to quantify their influence-may be included either in a single linear predictor or hierarchically.
Comprehensive spatial models should respect the complex potential state-wide correlation structure. Given the results in Sect. 3.1, it seems reasonable to focus on the river segments when building up a neighbourhood structure of the water supply stations based on their locations along the NRW river network.
Within our PerSpat project, a fixed correlation matrix of the water supply areas based on river dependency and the estimated station-to-area-supply proportions is initially used to apply an algorithm solving a realignment problem (Kohlenbach 2019) and also for the preliminary smoothing approaches. In the future, river-induced correlations may be estimated from the Ruhr data.
Spatial correlations caused by groundwater may be relevant, too, and, therefore, are possible extensions of the neighbourhood structure. The existence of several contamination sources should be taken into account, too.
With our data amount, all such models are reliable only for the river Ruhr, but generally applicable to all main rivers of the state. Our modelling approaches are transferable to other or more complex river systems and to other substances solute in water or to similar scenarios.
Predicted PFOA concentrations in drinking water will be used in further phases of our PerSpat project to estimate PFOA concentrations in plasma of NRW residents, depending on time and place of residence, and for the analysis of health related data from the state-wide birth registry.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.