Determination of groundwater threshold values: A methodological approach

Journal Pre-proof Determination of groundwater threshold values: A methodological approach Onur F. Bulut, Burcu Duru, Özgür Çakmak, Özgür Günhan, Filiz B. Dilek, Ulku Yetis PII:

S0959-6526(20)30048-2

DOI:

https://doi.org/10.1016/j.jclepro.2020.120001

Reference:

JCLP 120001

To appear in:

Journal of Cleaner Production

Received Date: 3 March 2019 Revised Date:

1 January 2020

Accepted Date: 3 January 2020

Please cite this article as: Bulut OF, Duru B, Çakmak Öü, Günhan Öü, Dilek FB, Yetis U, Determination of groundwater threshold values: A methodological approach, Journal of Cleaner Production (2020), doi: https://doi.org/10.1016/j.jclepro.2020.120001. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.

1

Determination of Groundwater Threshold Values: A Methodological Approach

2 Onur F. Buluta, Burcu Durub , Özgür Çakmakc, Özgür Günhanc, Filiz B. Dileka, Ulku

3

Yetisa*

4 a

5

Middle East Technical University, Department of Environmental Engineering, 06800

6

Ankara, Turkey b

7 8

c

Fugro Sial, Guvenevler Mahallesi, Farabi Sk. 40/4, 06690 Ankara, Turkey

Ministry of Agriculture and Forestry, General Directorate of Water Management, 06560

9

Ankara, Turkey

10

Abstract

11

This article is concerned with the establishment of natural background levels and threshold

12

values for naturally-occurring parameters in groundwaters in the absence of relevant and

13

accurate long-term spatial data. A new approach was developed and exemplified, adopting the

14

groundwaters of the Gediz River Basin, Turkey, as a case study. The available groundwater

15

monitoring data was from a one-year seasonal water quality monitoring campaign. The

16

approach used combines pre-selection, selection and statistical evaluation of the selected data

17

to eliminate outliers to determine natural background levels. The application of three different

18

statistical tools, namely, probability plot, 2-σ iteration, and distribution function methods,

19

resulted in different natural background level estimates. The 2-σ iteration method provided

20

the most conservative values for almost all the parameters. The use of this three-step

21

approach, which adopts different statistical methods, appeared to solve the limited data

22

availability challenges specific to groundwater contamination and improve natural

23

background level assessment and threshold value setting. Lessons learned from this study can

24

help policymakers to promote similar initiatives in other countries where groundwater quality

25

data is limited.

26

*

Corresponding author E-mail address: [email protected]

1

27

Keywords: Groundwater, pre-selection, outliers, threshold value, natural background level

28 29

Nomenclature

30

BRIDGE: Background cRiteria for the Identification of Groundwater thrEsholds

31

EC: European Commission

32

EQS: Environmental Quality Standards

33

EU: European Union

34

GWB: Groundwater Body

35

GWD: Groundwater Directive

36

LoQ: Limit of Quantification

37

MoFW: Ministry of Forestry and Water Affairs

38

NBL: Natural Background Level

39

REF: Reference Value

40

TBGW: Turkish bylaw on the Protection of Groundwater against Pollution and Deterioration

41

TV: Threshold Value

42

1. Introduction

43

The pressure on water resources that threatens water safety is a growing concern (Guo

44

and Bartram, 2019) due to population growth, industrialization, and agriculture (Ramírez et

45

al., 2017). In many parts of the world, the freshwater available is decreasing in quantity and

46

worsening in quality (Das et al., 2019) and groundwater reserves which are substantial

47

freshwater resources are diminishing (Niu et al., 2019) due to higher rate of water pumping of

48

aquifers than the natural recharge rate (Rodríguez-Huerta et al., 2019). Further pressure on

49

groundwaters is due to poor sanitation and weak infrastructure, as well as geogenic

50

contamination. Developing management strategies for the protection of groundwater and

51

assessing the risk they may cause for specific use are significant concerns worldwide as that

52

improves public access to safe drinking water ((Jenifer and Jha, 2018).

53

The European Groundwater Directive (2006/118/EC) on the protection of groundwater

54

against pollution and deterioration that is commonly referred to as the “Groundwater

55

Directive” (GWD) aims to establish specific measures to prevent and control groundwater

56

pollution. These measures include criteria for the assessment of good chemical status, the 2

57

identification and reversal of environmentally significant pollutant trends, and preventing or

58

limiting inputs of pollutants into groundwater. Accordingly, the GWD requires the use of the

59

groundwater quality standards defined in the Directive as well as threshold values (TVs) to be

60

set by following the application of the procedure and considerations set out in the Directive.

61

TVs are groundwater quality standards that are set by the EU Member States and several

62

other states to test whether groundwater bodies are at good status during the groundwater

63

body classification task. Article 3 of the GWD states that the Member States should use

64

prescribed groundwater quality standards for nitrates and pesticides, and locally derived TVs

65

for other pollutants that were identified as putting the groundwater body (GWB) at risk of

66

failing environmental objectives. The GWD provides a minimum list of pollutants that the

67

Member States are asked to consider when setting TVs. The Turkish bylaw on the Protection

68

of Groundwater against Pollution and Deterioration (TBGW) taking the GWD as the basis,

69

also requires the determination of TVs for the target chemical pollutants towards sustainable

70

management of groundwater bodies. Within this framework, the General Directorate of Water

71

Management carried out a pilot project for the development and implementation of

72

methodologies for the assessment of groundwater quantity and quality, in line with the above-

73

mentioned by-law (MoFW, 2017). Within the content of this project, all provisions of the by-

74

law were implemented, and a methodology was developed for setting TVs for the Gediz River

75

Basin.

76

The TVs need to be set so that the natural background levels (NBLs) of chemicals in

77

groundwater do not lead to poor status. Pollutants in groundwater are either caused by human

78

activities (Collin and Melloul, 2003) or by natural geochemical reactions between the

79

groundwater and the rock-forming the aquifer (Jenifer and Jha, 2018). NBLs should be based

80

on the monitoring data, which should be collected from monitoring points representing

81

different flow conditions and chemistry throughout the area and depth of the GWB. For

82

example, in confined aquifers, the natural chemistry of groundwater will be quite different

83

from that in unconfined aquifers. If there is limited groundwater quality data, it is suggested

84

that the NBL can rely on selecting a subgroup that does not show pollution and takes into

85

account groundwater flow and geochemical conditions in the GWB (GWD 2000/118/EC).

86

The essence of setting TVs lays in the determination of the NBLs, which are affected by

87

several factors; such as, chemical and biological processes in the unsaturated zone, the

88

residence time of water in the aquifer, recharge by rain, relations with the other aquifers and 3

89

water-rock interactions (Urresti-Estala et al., 2013). Although the first approach to NBL

90

estimation was based on a geochemical perspective in which the concentration data is

91

belonging to aquifers on which anthropogenic pressure does not exist, the typical approach

92

employed is by the statistical analysis of monitored water quality data (Edmunds and Shand,

93

2008). In cases where the concentration data is from an aquifer that is not pristine, the

94

geochemical approaches cannot certainly be used for the estimation of NBL. Instead,

95

statistical methods should be used to separate the natural population of geochemical data from

96

the anthropogenic one by assuming a normal or log-normal distribution (Molinari et al.,

97

2012). In line with this, the European research project entitled BRIDGE (Background cRiteria

98

for the Identification of Groundwater thrEsholds) suggests a methodology termed as pre-

99

selection, that is based on statistical evaluation of water quality data and the identification of

100

pristine groundwater samples as representative of NBL (Müller et al., 2006). The BRIDGE

101

methodology for setting NBL relies on two major steps. The data sets are first eliminated by

102

applying specific pre-selection criteria and then the 90th or 97.7th percentile of the pre-

103

selected data is chosen as NBL (Müller et al., 2006). As reviewed by Sellerino et al. (2019),

104

there are many other studies in the literature that combine pre-selection methods and

105

statistical NBL determination. For example, Parrone et al. (2019) applied the pre-selection of

106

uncontaminated samples by using nitrate/ammonia as the appropriate marker and adopted the

107

NBL setting based on a normal distribution. Molinari et al. (2012) applied component

108

separation and pre-selection methodologies and indicated that the estimated values of NBLs

109

by these two approaches are within the same order of magnitude.

110

Indeed, all the methods applied for the NBL setting rely on the elimination of the outliers

111

or extreme values assuming that the extreme values or outliers are not from the background

112

but another process or source, and thus the remaining data belong to the background. To this

113

end, several different statistical techniques have been used. For example, Preziosi et al. (2014)

114

used the probability plot method, Matschullat et al. (2000) used the 2-σ iteration technique,

115

Cidu et al. (2017) used the median absolute deviation method, Reimann et al. (2005) used the

116

box and whisker plot. As a one-step further approach, Sellerino et al. (2019) adopted both the

117

probability plot and pre-selection methods to assess NBL for As, F, Fe, and Mn and came up

118

with different NBLs, depending on the adopted method. They attributed this difference to the

119

variation in hydrogeochemistry in the groundwater body and reported that the NBLs

120

estimated for F and Mn by the probability plot method appear to be inappropriate, unlike for 4

121

As and Fe, probably due to the fact that all the concentration values are very high and the

122

geochemical anomaly is extended to the whole groundwater body. Therefore, the results from

123

the pre-selection method were considered as the final NBL. In line with this finding, Mollinari

124

et al. (2019) recently claimed that NBLs could have very different local values, due to the

125

occurrence

126

groundwater body, especially in large-scale reservoirs with an areal extent of several

127

thousands of square kilometers. Hence, they reported that the common practice of evaluating

128

the chemical status relying on a single NBL value could not be considered as realistic.

129

Accordingly, they proposed an approach for the estimation of local NBLs at

130

the borehole scale through the application of a geostatistically-based methodology to yield

131

exceedance probability maps. All these studies which rely on functional and comprehensive

132

datasets were satisfactory despite the use of a single statistical method for NBL determination

133

due to the availability of relevant and accurate long-term spatial data. However, developing

134

countries such as Turkey face significant challenges when setting TV, given the scarcity of

135

long-term groundwater quality data. Such long-term comprehensive groundwater quality data

136

is difficult to collect, mainly due to the high cost of monitoring campaigns.

of

diverse

petrographic provinces or redox

conditions within

the

same

137

An approach to overcome the data scarcity problem is to design techniques that can deal

138

with minimal data sets. The present study was directed by limited data availability and

139

developed a new approach to NBL determination for the TV setting. The methodology

140

developed was mainly based on the assessments used in the EU Member States (BRIDGE

141

Project) and the groundwater quality dataset available. As different than the BRIDGE

142

methodology which includes the steps of pre-selection and NBL setting with the use of

143

probability plot method, the NBL determination methodology developed in this study

144

included i) data pre-selection, ii) elimination of outlier data using box and whisker plots (data

145

selection), and iii) determination of NBL applying three different statistical techniques

146

(probability plot, 2-σ iteration and distribution function methods). More specifically, we used

147

the additional step of data selection and adapted a modified version of the statistical NBL

148

setting with the use of three different statistical techniques to produce more realistic and more

149

conservative NBLs. Data selection was added to the methodology; because the pre-selection

150

criteria adopted in the BRIDGE Project were not fully applicable in the present study due to

151

the limited number of data available and therefore, there was a need for better elimination of

152

outliers. NBL setting using three different statistical techniques was considered because the 5

153

distribution of data was not always normal. Also included is a comparison among the NBLs

154

derived by the different statistical techniques to evaluate how far the NBLs estimated are

155

sensitive to the chosen statistical technique. By combining the pre-selection method with the

156

above-mentioned statistical methods and by verifying its applicability at the site, it is aimed to

157

have a conservative and flexible approach that can be adapted according to the availability of

158

data. After NBLs are determined, reference (REF) values are chosen and finally, TVs are

159

determined by comparing NBLs and REFs.

160

2. Materials and Methods

161

In this section, a description of the study area, data set available, pollutants considered,

162

and NBL determination performed during the TV determination for the contaminants of

163

concern are provided.

164

2.1. Study Area

165

The study area is the Gediz River Basin (Figure 1), which is among the nine high priority

166

river basins in Turkey, as identified by the Action Plan on Groundwater Management in 2013.

167

The Basin drained by the Gediz River has an approximate area of 17,500 km² and hosts very

168

fertile agricultural lands, animal husbandry activities, organized industrial sites, high potential

169

geothermal fields, variety of mineral deposits, in addition to the densely populated

170

settlements. The Gediz River has a length of 400 km, with an average discharge of 60 m3/s at

171

its mouth. About two-thirds of the basin area stay natural, free from anthropogenic activities,

172

mostly located in the northern and northeast parts of the basin (MoFW, 2017). The other part,

173

which is the downstream area of the Gediz River basin, is home to a variety of industrial

174

activities such as textile, food processing, leather, dairy, meat, poultry processing, and

175

agricultural vehicle manufacturing. The discharges from such plants significantly raise the

176

levels of heavy metals and organic micro-pollutants in the river network. The Gediz River

177

Basin grows a variety of agricultural products that mainly include grape, olive, cherry,

178

tomato, walnut, and cotton.

179

Geology of the Gediz Basin was studied in detail in a previous study by the State

180

Hydraulic Works of Turkey and summarized in the Hydrogeological Investigation Report for

181

the Gediz Basin (SHW, 2015). According to this report, basement rocks of the basin are made 6

182

up of metamorphic rocks of Menderes Massive. Paleozoic units are unconformably overlain

183

by Mesozoic schists intercalated with metaconglomerates. These units then overlap with

184

carbonates, while the transition zone is identified by alternating dolomite, quartzite, and

185

calcschists. Massive dolomites overlying these alternating units are located beneath very thick

186

(reaching to 1500 m) massive marbles. İzmir-Ankara zone, made up of ophiolite and flysch

187

units, is situated on top of Menderes metamorphics by the thrust faults. These units are

188

unconformably overlain by terrestrial and lacustrine sequences of Neogene sedimentary units,

189

volcanic and igneous units. Quaternary basalts and alluvium units cover all the above-

190

mentioned base units of the Gediz Basin. Basalts in the Kula region are well-known for about

191

80 volcanic cones of lava and tephra. Other Quaternary units of the basin are the uncemented

192

alluvium, talus, fan and terrace deposits.

193

Moreover, hydrogeology of the Gediz Basin was also studied within the scope of that

194

report, where hydrogeological properties of the geological units and their corresponding

195

groundwater-abstraction potential (based on their specific discharge, hydraulic conductivity,

196

transmissivity, well and spring yields, etc.) were determined. Based on the results of these

197

studies, aquifers of the basin were identified. The karstic rock groundwater bodies made up of

198

marble, limestone, dolomite, and travertine are described as aquifers providing significant

199

amounts of groundwater. In addition to these karstic units, granular units, which are

200

commonly observed in the basin and are deposited as alluvial sediments, alluvial fans, cones,

201

and slope debris; also have significant groundwater potential. In the basin, groundwater can

202

be obtained with high rates from Neogene aged clastic rock mass, depending on the sandstone

203

and conglomerate levels it contains. Similarly, Neogene aged volcanic rocks in the basin can

204

provide groundwater at the regional and local scale where they have secondary porosity. In

205

addition to these units, clayey limestone units of Neogene limestones in the basin are also

206

used to supply groundwater locally. However, in the regions where the clay content is high,

207

specific capacities of wells drilled in this unit are reduced. On the other hand, Paleozoic

208

metamorphic rocks and Mesozoic flysch units, having very low specific capacities, were

209

defined as the units not having the potential to be classified as aquifers. Furthermore, in this

210

former study, all the aquifers were grouped according to their groundwater-abstraction

211

potentials, as follows: Neogene clastic rocks, volcanic rocks, and clayey limestones of the

212

basin were classified as lower-yield aquifers of limited groundwater potential; while

213

Paleozoic marbles, Mesozoic and Neogene limestones and uncemented units (alluvium, 7

214

alluvial fans, talus) were classified as higher-yield aquifers of significant groundwater

215

potential.

216

217 218

Figure 1. Location map of the Gediz River Basin (Karaaslan, 2017)

219 220

8

221

2.2. Available Dataset

222

Previous work on water quality monitoring in the basin was limited to parameters such as

223

Ca, Mg, Na, electrical conductivity, etc. Within the framework of the present study (MoFW,

224

2017), for the first time, a monitoring program including many metals, metalloids, anions was

225

developed for the basin to determine the quality of groundwater bodies. Considering the

226

duration and the scope of this project, it was deemed appropriate to monitor three periods

227

during the project in such a way to represent both wet conditions and dry conditions. During

228

the design of the sampling network, the locations were selected in such a way that they

229

represent the quality of groundwater bodies, considering the locations of all wells and springs

230

in the basin.

231

The sampling network included wells and springs with variable spacing, covering the

232

entire basin. Some of the wells/springs were located far from the anthropogenic pollutant

233

sources to obtain some pristine water samples, while some others were located nearby

234

anthropogenic pollutant sources, as the monitoring activity was run not only for NBL

235

assessment but also overall groundwater quality assessment. Also, groundwater flow

236

directions were taken into consideration, and possible areas that may be impacted by these

237

anthropogenic pollutant sources were considered to reveal the effects on groundwater quality.

238

A total of 110 sampling locations representing 71 out of 76 groundwater bodies within the

239

basin were identified and samples were collected (Figure 2). There were no wells/springs

240

suitable for sampling in five of the groundwater bodies, and no chance of drilling new wells

241

within the scope of the project.

242

2.3. Pollutants for NBL derivation

243

According to the Turkish Legislation, for the groundwater bodies which are classified to

244

be at risk, it is required to set TVs at the most appropriate scale (national, river basin district

245

or groundwater body) for each parameter that causes the groundwater body to be classified as

246

at risk. In Turkey, for the common groundwater pollutants of nitrates and pesticides, there are

247

quality standards set at the national scale by TBGW. Therefore, within the scope of the

248

present study, only 29 pollutants (excluding nitrates), which are not purely anthropogenic and

249

that may cause groundwater bodies in the Gediz Basin to be classified at risk, were taken into

250

consideration. In doing this, the data collected for a total of 29 pollutants were firstly 9

251

processed to estimate the average concentration for each parameter. For the calculation of

252

average concentrations, values below the limit of quantification (LoQ) were set to half of the

253

value of the LoQ concerned as suggested by the Guidance Document No. 19 of the European

254

Commission (EC, 2009). The average concentrations were then compared with the

255

corresponding LoQs, to produce an average dataset that represents the average groundwater

256

characteristics. In case the calculated mean value for a parameter exceeds its LoQ at least at

257

three sampling stations, this parameter was considered as posing a risk; therefore, its TV has

258

to be determined. There were seven parameters (F-, CN-, Ba, Be, Sb, Ti, and Ag) which were

259

classified as not posing a risk. As a result, a list of 22 naturally occurring parameters given in

260

Table 1 was formed to take into account in assessing NBLs and in setting TVs.

261

Supplementary Materials - Part I presents the monitoring data collected for 22 posing risk and

262

7 non-risk posing parameters separately.

263 264

Table 1. Parameters for which NBL and TV are determined Group

Substances

Ions

Cl-1, SO4-2, S-2, PO4-3-P

Metals

Cd, Hg, Cu, Zn, Fe, Co, Mn, Mo, Ni, V, Cr, Pb, Na, Al

Metalloid s Other

As, B, Se Electrical conductivity

10

265 266

Figure 2. Position of the Quality Monitoring Points in the Study Basin (MoFW, 2017)

267 268

2.4. Determination of Threshold Value

269

As also stated earlier, when dealing with large scale aquifer systems having data scarcity,

270

it could be convenient to analyze available data monitored through statistical analysis

271

methods. In this respect, the method presented in Figure 3 was followed in developing TVs.

272

This method is mainly composed of i) NBL assessment which includes data pre-selection,

273

outlier elimination (or data selection) by the box and whisker plots, and NBL estimation using

274

three different statistical methods, ii) determination of REF value, iii) TV setting based on the

275

comparison of the NBL and the selected corresponding REF value.

276

2.4.1. Assessment of Natural Background Level

277

Pre-selection: Within the framework of the BRIDGE project, Wendland et al. (2005,

278

2008) proposed a pre-selection as a simplified approach for NBL determination. Pre-selection

279

is applied where a limited set of quality data is available and therefore, it is not possible to 11

280

adopt a component separation method. In the component separation method, the observed

281

concentration frequency distribution is fitted by the superimposition of two individual

282

distributions that represent the natural and the influenced component. Once the distribution of

283

the natural component is assessed, the related data are used to estimate the NBL (Müller et al.,

284

2006).

285

The pre-selection methodology involves selecting sampling locations that meet specific

286

criteria for the exclusion of samples affected by human influence. The adopted criteria in the

287

present project, encompassed the following constraints: (a) locations where chloride

288

concentration >1000 mg/L, (b) locations where nitrate concentration >50 mg/L, and (c)

289

locations where there is geothermal pressure are assumed to be contaminated and removed

290

from the data set. The other exclusion criteria that would ideally be adopted (e.g., redox

291

conditions, anaerobic conditions) cannot be applied, since it would be more appropriate not to

292

apply, given that there are minimal numbers of available sampling points that could be used.

293

Likewise, the nitrate concentration criterion of <10 mg/L suggested by the BRIDGE project

294

was found to be inapplicable, as the Gediz Basin is with intensive agricultural activities and

295

most of the monitoring stations have high average nitrate concentrations.

296

Selection - Box and Whisker Plots: After pre-selection, distribution of the data for each

297

parameter was examined using box and whisker plots for determining whether the distribution

298

is skewed and whether there are potential unusual values or outliers in the data set. This step

299

forms the selection of the data to be used in the determination of NBL using statistical

300

methods. The values in the data set under LoQ values were assumed to be equal to LoQ/2 and

301

included in creating box and whisker plots. Box and whisker plots are used to visualize the

302

distribution of the dataset, and to identify outlier concentrations that are likely to be due to

303

man-made impacts. The median is the centerline of the box plot. The upper limit of the box

304

represents the 25th and the lower limit of the box represents 75th percentile values, while the

305

whiskers or the upper and lower bars represent the outer limits of the dataset (e.g., plus or

306

minus two standard deviations from the mean value). The maximum value of the data set

307

should not be more than 1.5 times the difference between 75% and 25%. With this principle,

308

outliers are determined.

309

NBL Setting: Subsequently, the NBL is estimated based on the modified distribution of

310

the concentration data, applying three different methods, which are probability plot, 2-σ

12

311

Iteration and distribution function methods. After the determination of NBL by the

312

application of these methods, the lowest value was considered as the final NBL.

313

314 315 316

Figure 3. Method Used for Setting Threshold Values (EQS: environmental quality

317

standard for surface waters, NBL: natural background level, REF: reference value, TV:

318

threshold value, WHO: World Health Organization)

319 320

The Probability plot method is also proposed by the BRIDGE project (Müller et al.,

321

2006). After pre-selection and selection are made, the probability plot for the data remaining

322

is prepared (using MATLAB software) and the value of covering 90% or 97.7% of the dataset

323

is evaluated and accepted as NBL. The selection of the percentile depends on the size of the

324

data set. For small data sets (<60 data), 90% is used, while for large data sets, 97.7% is used

325

(Hinsby et al., 2008). 13

326

The 2-σ iteration method treats the data set by iteration until a normal distribution is

327

constructed around the mode of the original data. This technique calculates the mean and the

328

standard deviation of this normal distribution and establishes a range, of which the values

329

outside are eliminated. Then, the Lilliefors test is applied to check the normal distribution of

330

data. During the Lilliefors test, the t-test is performed and the calculated t-value is compared

331

to the t-critical value (α value is 0.05 or 95% confidence level). If the calculated t-value is

332

smaller than the t-critical value, it is decided that the obtained distribution is normal. In this

333

case, this method is deemed suitable for calculating the NBL (Urresti-Estala et al., 2013), and

334

the mean+ 2-σ value is taken as NBL. Nakic et al. (2007) pointed out that this technique is

335

appropriate for calculating TVs as the upper limit of the NBL. The suitability of this method

336

for NBL assessment depends on the frequency distribution of the data and the properties of

337

the parameter (Urresti-Estala et al., 2013). In cases where the difference between the

338

minimum and maximum concentration of a substance is very high, this method is indicated to

339

be unsuitable.

340

The Distribution function method, which is based on the concept that anthropogenic

341

factors generate an enrichment of the parameter values, works similarly as the 2-σ Iteration

342

method and looks for the normal distribution. The difference is that the concentration data

343

above the median are removed from the data set and the concentration data between the

344

minimum and median are considered to be free from human effects (Matschullat et al., 2000)

345

and therefore represent NBLs (Urresti-Estala et al., 2013). The Lilliefors test is applied to the

346

data between the minimum and the median values to check for normal distribution. The

347

resulting mean+ 2-σ value is taken as NBL. Like for the 2-σ Iteration method, the suitability

348

of this method depends on the frequency distribution of the data and the nature of the

349

parameter.

350

2.4.2. Determination of Reference Value

351

As shown in Figure 3, following the determination of NBLs, it is necessary to determine

352

REF values to be considered for establishing the TVs. Generally, when water quality needs of

353

groundwater-dependent ecosystems are not known, surface water mean annual environmental

354

quality standards (EQS) or drinking water standards are used as proxies for REF values

355

(Danielopol et al., 2003; Hose, 2005). Considering that groundwater bodies in the Gediz

14

356

Basin are used as drinking water and/or irrigation water, REF values have been set through

357

the following steps:

358

•

If groundwater is used for more than one purpose, i.e., for drinking water and

359

irrigation water then the most stringent standard is used as REF. If only one of these

360

standards is available, the available value is considered to be the criterion value.

361

•

In the absence of both drinking water and irrigation water criteria specified in

362

regulations, the relevant EQS value set out for surface waters (if any) is considered to

363

be the REF value.

364

•

For the parameters that do not have drinking water/irrigation water criteria and EQS

365

value specified in the regulations, according to EU Technical Guidance Document No.

366

27, EQS value can be derived or World Health Organization (WHO) drinking water

367

criteria can be used.

368

2.4.3. Setting Threshold Value

369

In the process of establishing TVs, like for the calculation of the NBL, all alternative

370

approaches applied in EU countries and also the one by the BRIDGE project (Müller et al.,

371

2006) were examined, and the most appropriate methodology for the conditions existing in

372

our country was determined as follows:

373

If NBL / REF ≥ 1; TV = NBL

374

If NBL / REF<1; TV=REF

375

3. Results

376

This section presents groundwater chemistry of the study area and the NBL and TV

377

setting results obtained with the application of the suggested methodology for the pollutants

378

in the Gediz River Basin under study.

379 380

3.1. Groundwater Chemistry of the Study Area

381

Natural groundwater chemistry is influenced by several factors such as aquifer lithology

382

and geochemistry, groundwater flow paths, groundwater residence time, and rainfall 15

383

chemistry. When it rains, the rainwater infiltrates through the soil, subsoils, and bedrock, and

384

eventually reaches the groundwater. The degree to which groundwater becomes mineralized

385

depends on several factors: the aquifer mineralogy, the permeability, and nature of

386

groundwater flow, the presence and nature of overlying deposits, the pH and redox

387

conditions, groundwater flow paths and the groundwater residence times (Tedd et al., 2017).

388

Hydrogeochemistry of the basin was analyzed based on the results from the monitoring

389

program forming the basis of this research (MoFW, 2017). Owing to its complex geology, the

390

water chemistry of the Gediz River Basin is also very diverse. Piper diagrams were used to

391

estimate the geochemical characteristics of the water samples. The Piper plot representing all

392

the groundwater samples indicates that the overall geochemical composition of groundwater

393

in the Gediz Basin (Figure 4) falls into calcium bicarbonate type reflecting the dominant

394

geology of the aquifers in the basin. On the other hand, individual units are highly diversified

395

in terms of their hydrogeochemistry, even the samples were taken from the same geological

396

unit, for instance, alluvium, display a wide diversification (Mg-Ca-HCO3, Ca-Na-Mg-HCO3,

397

Na-Ca-Mg-HCO3, Ca-Na-HCO3 ve Ca-Mg-HCO3), depending on the water-rock interactions

398

and anthropogenic impacts. Moreover, it should also be noted that the Gediz Basin is very

399

rich in geothermal resources and mineralization, which also results in the natural enrichment

400

of several parameters in groundwater samples, such as Fe, Mn, As, etc.

16

401 402

Figure 4. Piper plot of groundwater in the Gediz River Basin

403 404

3.2. Natural Background Level Setting

405

In assessing NBLs, the first step was data pre-selection. When the pre-selection criteria

406

mentioned in Section 2.4.1 were applied to the original data set presented in the

407

Supplementary Materials - Part I, a subset of the original data, which is classified as

408

preselected data was formed for 22 parameters. As shown, almost half of the original data was

409

eliminated.

410

At the data selection step, data after pre-selection was fitted to a normal distribution on

411

the assumption that the range of NBL data should follow this distribution once all polluted

412

samples were removed as outliers. With the elimination of outliers that are likely to be due to

413

man-made impacts using box and whisker plots, the distribution of data was understood

414

better, and the data was made ready for further statistical analysis and identification of the

415

NBLs. Supplementary Materials - Part II presents box and whisker plots for only 19 of the

416

considered 22 parameters, as there was no data left after pre-selection for S-2, B, and Se. 17

417

As can be seen from Table 2, there was a remarkable decrease in the number of sampling

418

points that provided acceptable data for NBL assessment, and in turn, in the dataset after the

419

selection using the box and whisker plots. The number of sampling points that can be

420

considered in the NBL assessment was about half of the sampling points from which water

421

quality data were collected during the monitoring campaign. Table 2 also indicates that there

422

was a serious change in the statistical characteristics of the water quality data following the

423

data pre-selection. For almost all the parameters, 75th percentile values were found to

424

decrease as compared to the original data. These results have confirmed that the groundwaters

425

in the Gediz River Basin are under a very serious anthropogenic stress, and therefore the task

426

of NBL setting in the presence of the data scarcity problem is a challenging task for this river

427

basin.

428

After removing any identified outliers with the use of box and whisker plots, the

429

remaining data were used for the evaluation of NBLs for each of the parameters, applying

430

probability plot, 2-σ iteration, and distribution function methods to remove any further

431

outliers. The NBL for each method was based on a percentile of the dataset (Table 3). The

432

percentiles selected for the probability plot method are based on the values suggested by the

433

BRIDGE project (Hinsby, 2008). For the 2-σ iteration and the distribution function methods,

434

the 95th percentile of the dataset was selected based on the method used by Urresti-Estala et

435

al. (2013).

18

436

Table 2. Number of groundwater samples in sample points and statistical parameters for 22 pollutants before and after pre-selection Original Data Parameter

Unit

Number of Sampling Points

Min

Max

25th

Ions Cl-1 SO4-2 S-2 PO4-3-P

mg/L mg/L mg/L µg/L

110 110 11 108

5.48 2.01 0.01 0.50

100.57 550.95 0.60 1521.57

12.37 15.81 0.01 5.64

Cd Hg Cu Zn Fe Co Mn Mo Ni V Cr Pb

µg/L µg/L µg/L µg/L µg/L µg/L µg/L µg/L µg/L µg/L µg/L µg/L

81 54 102 96 97 100 101 79 95 87 100 83

0.03 0.03 0.50 1.25 15.00 0.07 0.50 0.50 0.83 0.50 0.50 0.50

0.14 3.64 81.08 456.56 7165.50 7.14 1420.32 50.68 70.31 46.03 33.31 23.25

0.03 0.03 1.20 2.44 33.25 0.10 1.80 0.50 1.14 0.72 1.10 0.69

Na As Al B Se

µg/L µg/L µg/L µg/L µg/L

110 74 110 15 11

350.8 2.50 1.48 250.00 2.50

618418.7 3252.50 3481.08 7176.47 15.83

3058.2 2.50 2.65 250.00 2.50

110

33.15

1523.33

518.65

Electrical conductivity

µS/cm

After Preselection Percentile

Number of Sampling Points

Min

Max

19.03 37.29 37.70 75.78 0.01 0.01 19.10 47.38 Metals 0.06 0.07 0.03 0.06 2.01 3.41 5.74 13.57 90.76 533.90 0.17 0.29 6.69 61.48 1.52 3.12 1.53 2.71 2.22 6.26 1.89 4.81 1.57 2.59 Metalloids 5943.6 11205.4 9.60 35.73 6.56 18.30 250.00 250.00 2.50 2.50 Other

62 62 3 61

5.48 5.47 0.01 0.50

39 28 57 50 53 55 57 42 51 44 55 49

650.00

50th

437 19

75th

839.42

Percentile 25th

50th

75th

67.72 394.71 0.43 1521.57

11.20 12.03 0.01 3.77

14.80 28.22 0.01 14.83

24.89 60.32 0.01 32.87

0.03 0.03 0.50 1.25 15.00 0.07 0.50 0.50 0.83 0.50 0.50 0.50

0.13 1.20 81.08 285.36 7165.50 3.91 1049.60 50.68 20.43 46.03 29.09 23.25

0.03 0.03 0.95 2.38 27.36 0.10 1.70 0.50 1.06 0.50 0.90 0.87

0.05 0.03 1.85 5.17 96.48 0.15 8.23 1.13 1.36 1.96 1.68 1.61

0.06 0.06 3.05 12.65 533.90 0.25 35.63 2.82 2.47 6.73 4.19 2.47

62 40 62 5 5

350.8 2.50 1.48 250.00 2.50

55438.6 1348.16 2856.05 3882.33 6.07

2121.8 2.50 2.48 250.00 2.50

4532.5 9.53 5.43 250.00 2.50

8058.5 33.81 13.81 250.00 2.50

62

223.97

1202.67

452.76

573.42

782.33

438

The probability plots for all 19 parameters are given in the Supplementary Materials –

439

Part III. A straight line in a probability plat indicates that the data fit the normal distribution

440

and a deviation from this straight line could indicate an upper limit for the background

441

population, hence a possible threshold between natural and influenced concentrations

442

(Preziosi et al., 2014). A review of the probability plots indicated that most of the parameters

443

follow near to normal distribution except Fe, Mn, V, and Hg, which tend to exhibit relatively

444

a poorer distribution. The probability plots for Al, As, Cr, Zn, Cd, Co, Mo, and Ni show a

445

bimodal distribution. Among those with a high frequency of results for LoQ concentrations

446

are As, V, Cd, and Mo. The dataset for these four parameters includes a high number of

447

results reported at LoQ/2 values. The inclusion of a high number of LoQ results is likely to

448

skew the selection of the 90th and 97.7th percentiles away from actual NBLs when using

449

these probability plots. As presented in Table 4, NBLs were calculated for all 19 parameters

450

which can pass the data selection step.

451 452

Table 3. Percentiles used with different statistical methods Method

Percentile used for NBL

Probability plot < 60 data points

90

Probability plot > 60 data points

97.7

2-σ iteration method

95

Distribution function method

95

453 454

Following the evaluation of NBLs by probability plots, the 2-σ iteration method was

455

applied. These two methods were found to be applicable for the parameters having sufficient

456

data that fit into a normal distribution. As indicated in Table 4 and Supplementary Materials-

457

Part III, for the parameters with relatively insufficient data (Cd, Hg, Fe, Mo, V, As), this

458

method could not be applied. For these parameters, the Lilliefors test indicated that the

459

method did not work. It was due to the occurrence of a high number of results reported at

460

LoQ/2 values (for As, V, Cd, Mo) or a poorer distribution in probability plots (for Fe, V, Hg).

461

Similarly, for 3 of the 6 parameters mentioned above, namely As, Hg and Cd, the

462

distribution function method could also not be applied as indicated by the Lilliefors test

20

463

results. The distribution function plots for 16 parameters are given Supplementary Materials-

464

Part III.

465 466

Figure 5 presents the results from the application of these three statistical methods to

467

electrical conductivity data from the basin as an example. As can be depicted from Table 2,

468

the electrical conductivity data points sitting outside the 223.97 to 1203.67 µS/cm range were

469

considered as outliers after the pre-selection step. The remaining data had 25, 50 and 75th

470

percentile values of 452.76, 573.42 and 782.33 µS/cm, respectively. When further outlier

471

elimination with the use of the box and whisker method was applied, it appeared that there are

472

no outliers to eliminate in the data. Figure 5 presents the probability plot prepared for

473

electrical conductivity. As shown, the distribution of the electrical conductivity data was

474

almost normal with a 97.7th percentile value of 1186 µS/cm. When the 2-σ iteration method

475

was applied to the data after selection, its suitability was confirmed with the Lilliefors test

476

result, giving the t value of 0.042 and the t-critical value of 0.122. This method yielded an

477

NBL value of 882.1 µS/cm. Similarly, the distribution function method was proven to be as

478

suitable as it was evident from the Lilliefors test results with a t value as 0.024 and t-critical

479

value of 0.114. The resulting NBL value was 881.9 µS/cm. Thus, three methods applied

480

provided similar NBLs and the lowest value of 881.9 µS/cm was selected as the NBL for

481

electrical conductivity, considering that choosing the lowest is a more conservative approach.

482 483

Table 4. NBL, REF, and TV of the parameters

Distribution function

NBL Final

Unit

2-σ İteration

Parameters

Probability plot

NBL REF

32.1 0.91 78.9 NA 48.7

20.91 0.59 42.5 NA 24.0

22 0.62 56.4 NA 26.2

20.9 0.59 42.5 NA 24.0

106.52 3 2504 0.0053 2003

TV

106.5 3 250 0.005 200

Ions Cl-1 SO4-2 S-2 PO4-3

mg/L mEq/L mg/L mg/L µg/L

21

Distribution function

NBL Final

Metals Cd Hg Cu Zn Fe Co Mn Mo Ni V Cr Pb Metalloids Na

2-σ İteration

Parameters

Probability plot

NBL REF

µg/L µg/L µg/L µg/L µg/L µg/L µg/L µg/L µg/L µg/L µg/L µg/L

0.07 0.06 4.8 13.0 914.7 0.26 57.6 3.7 2.6 10.4 5.3 3.0

3.2 3.6 0.2 4.0 1.6 2.2 3.1

2.8 7.0 170.1 0.2 8.4 1.7 2.1 4.1 3.2 3.2

0.07 0.06 2.8 3.6 170.1 0.2 4.0 1.7 1.6 4.1 2.2 3.0

54 14 2002 20002 2004 502 504 102 204 1002 504 104

5 1 200 2000 200 50 50 10 20 100 50 10

µg/L

8518.7

7581.4

7581. 4 39.7 3.6 NA NA

69000

69000

104 2004 10004 104

39.7 (53)5 200 1000 10

881.9

7002

882

Unit

As µg/L Al µg/L B µg/L Se µg/L Other Electrical conductivity µS/cm 1

484

39.7 18.5 NA NA

3.6 NA NA

8192. 9 8.7 NA NA

1186

882.1

881.9

NBL values considered as final NBL are in bold,

2

TV

2

Irrigation water standards,

485

3

Environmental quality standards, 4Drinking water standards, 5TV for As was set by the Ministry as

486

equal to the environmental quality standard, NA: not applicable due to no data left after pre-

487

selection.

22

488

489 490

Figure 5. (a) 2-σ Iteration (b) Probability (c) Distribution function plots for electrical conductivity 23

491

3.1. Threshold Value Setting

492

TVs determined for the parameters considered are presented in Table 4. As presented in

493

this table, REF values considered in deriving TVs are mostly drinking water standards as the

494

relevant authority aims to protect groundwaters in the basin at this level. For two of the

495

parameters (S-2 and PO4-3-P) EQS were considered, as there are no drinking water standards

496

or irrigation water standards for these ions.

497

According to the methodology applied for calculating final TVs, specified in the last

498

column of this table, TVs obtained using EQS as REF are mathematically rounded to a lower

499

integer. When NBL is greater than REF value, it is mathematically rounded to an upper

500

integer.

501

4. Discussion

502

Table 4 presents NBLs that were evaluated for the parameters considered. As can be seen,

503

in most cases, different methods yielded different NBLs. The differences between the

504

estimated NBLs are not mostly negligible. The variation was one order of magnitude for Al,

505

Mn, Zn, and V, which are with skewed distributions. As can be depicted from Table 4, this

506

difference was also more pronounced for the parameter of Fe, possibly due to a poorer fit of

507

the data into a normal distribution. What is also seen from Table 4 is that the methods of the

508

2-σ iteration and distribution function are in better agreement (except for Na), whereas the

509

probability plot method provided quite distinctly higher NBLs for almost all the parameters.

510

Then, it could be stated that for TV setting, the probability plot method was the least

511

conservative one and the 2-σ iteration and distribution function methods were more

512

conservative. Indeed, this could be because, in the case of the probability plot method, data is

513

not processed until a normal distribution is obtained, and therefore it is not satisfactory in

514

assessing NBL for datasets that are skewed and/or away from the perfect normal distribution.

515

The other two methods, however, process data until a normal distribution is obtained. So, it

516

seems that they yielded more reliable NBLs, in general. These observations were following

517

Uresti-Estala et al (2013) who stated that the NBLs obtained for the carbonate aquifers with

518

the 2-σ iteration and distribution function methods produce very comparable results because

519

the dominant geology of the Gediz Basin is carbonate (Figure 4).

24

520

On the other side, the applicability of the 2-σ iteration and distribution function methods

521

was limited by the skewness of data and also by the dispersion of data over a wide range of

522

concentration values. For example, in the case of As, which is a significant problem in the

523

Gediz Basin due to its natural presence at high levels in the groundwater (Sec 2.1.2), NBL for

524

As could not be determined with these two methods. This finding is not in agreement with

525

Urresti-Estala et al. (2013), who indicated that the 2-σ iteration method is more suitable when

526

the nature of the aquifer is responsible for high concentrations of parameters. The reason for

527

this contradiction was thought to originate from the skewed distribution of As data due to the

528

occurrence of high number of results below LoQ, as well as due to very wide range of As data

529

(
530

critical value. Urresti-Estala et al. (2013) also stated that the 2-σ iteration method is not

531

appropriate when the data set is dispersed and there are a large number of smaller values in

532

the data set, confirming the above-mentioned attribution. As also stated by them, the

533

suitability of these methods varies depending on the frequency and distribution of data. In

534

cases where the values are very variable, these methods are not suitable because the median or

535

mode of data is high and the number of data eliminated is high. In another study, supporting

536

Mapoma et al. (2016) could not implement these two methods for the parameter of As,

537

claiming that hydrochemistry of groundwater has a significant effect on As mobilization and

538

there are high numbers of As concentrations below LoQ in the As data set. Similarly, the 2-σ

539

iteration method and the distribution function method could not be implemented for Cd and

540

Hg because the calculated t-values were higher than the t-critical values.

541

On the other hand, for some parameters (Fe, V, and Mo), the 2-σ iteration method was not

542

found suitable as the Lilliefors test could not be passed while the distribution function method

543

was implemented successfully. For these parameters, many of the values in the dataset were

544

the same (i.e., LoQ/2) and the range of data was very wide (e.g., for Fe; min:
545

µg/L, max: 7165.5 µg/L). So, the above assessment for As is also valid for Fe. However,

546

unlike for As, although the 2-σ iteration method cannot be applied, the distribution function

547

method passed the Lilliefors test. It could be due to the large numbers of smaller values in the

548

dataset of Fe, V, and Mo as compared to As. Another reason could be the fact that the

549

presence of these metals is governed by different processes in the groundwater. The presence

550

of As can be related to water-rock interaction (Nath et al., 2018) with the volcanic formations

551

(Vivona et al., 2007) and/or to localized geothermal fluids (Angelona et al., 2009). In contrast, 25

552

the high values, for example, of iron may be related to the dissolution of Fe oxides/hydroxides

553

from the underlying alluvial sediments, as well as to its enhanced mobility by the anthropic

554

activities (Preziosi et al., 2014). On the other hand, mobility of Fe is naturally limited by the

555

strong buffering capacity of carbonates and by high cation exchange capacity of the organic

556

matters (Papadopoulou-Vrynioti et al., 2014).

557

Also, a particular emphasis should be given by the policymakers to the differences in

558

percentiles used with the different statistical methods. The difference between 90th, 95th and

559

97.7th percentiles can be very high for positively skewed distributions such as those of As, Fe

560

and Mn. In general, the different percentiles are proposed by different methods as per the

561

degree of uncertainty of the dataset used and the degree of knowledge of the

562

hydrogeochemical system and human influences. For example, Gemitzi (2012) used the

563

97.7th percentile during setting NBL in the plain area, and the 90th percentile in the coastal

564

area, where the data remained more limited after preselection. So, the elimination of specific

565

percentiles, such as 5%, 10% of the upper tail of the data set, needs to be justified depending

566

on the conditions and the practical effects of the degree of uncertainty of the dataset used in

567

the derivation of legally binding TVs should be carefully assessed.

568

5. Conclusions

569

This study presented a methodology to enhance the previously proposed threshold value

570

setting approach to deal with the data scarcity issues prevalent in developing countries. When

571

the methodology was implemented for the Gediz River Basin, it was seen that the

572

implementation of the three-step method developed that included pre-selection, selection, and

573

natural background level setting using probability plot, 2-σ iteration, and distribution function

574

methods provided more accurate natural background level estimates.

575

The pre-selection method proposed by prior studies was not applicable as it left no data to

576

process in the absence of relevant and accurate long-term spatial data. Applying the pre-

577

selection with a restricted set of criteria followed by data selection to eliminate outliers

578

statistically was satisfactory. It was seen that data selection by this approach before natural

579

background level assessment strengthened the process.

580

It was also concluded that the choice of statistical method to use in natural background

581

level assessment is a critical issue, as the results from different statistical methods may be 26

582

very different. The geological nature of aquifer, nature of pollutants, and the interaction

583

between them influence the suitability of the statistical method to use. Therefore, when there

584

is limited data, the integration of alternative statistical methods allowing the selection of more

585

confident lowest natural background level appears as a more robust and confident approach.

586 587

Acknowledgment

588

This study was performed within the scope of the project titled “Developing and

589

Implementation of Methodologies/Methods for the Determination and Assessment of

590

Groundwater Quantity and Quality: Gediz Basin Pilot Study” conducted by Fugro-Sial

591

Geosciences Consulting and Engineering under the authorization of the Ministry of

592

Agriculture and Forestry. The authors thank the General Directorate of Water Management

593

for their funding and guidance during the study.

594 595

References

596 597 598

Angelone, M., Cremisini, C., Piscopo, V., Proposito, M., Spaziani, F. 2009. Influence of hydrostratigraphy and structural setting on the arsenic occurrence in groundwater of the Cimino-Vico volcanic area (central Italy). Hydrogeol. J. 17, 901-914.

599 600 601

Cidu, R., Biddau, R., Lorrai, M., Mulas, M. G. 2017. Assessing Background Values of Regulated Parameters in Groundwater Bodies of Sardinia (Italy). Procedia Earth Planet. Sci. 17, 205-208.

602 603

Collin, M. L., Melloul, A.J. 2003. Assessing groundwater vulnerability to pollution to promote sustainable urban and rural development. J. Clean Prod. 11(7), 727-736.

604 605

Danielopol, D.L., Griebler, C., Gunatilaka, A., Notenboom, J. 2003. Present state and future prospects for groundwater ecosystems, Environ. Conserv. 30(2), 104-130.

606 607 608

Das, R., Laishram, B., Jawed, M., 2019. Perception of groundwater quality and health effects on willingness to procure: The case of upcoming water supply scheme in Guwahati, India, J. Clean Prod. 226, Pages 615-627

609 610 611 612 613 614

EC, 2009. Common Implementation Strategy for the Water Framework Directive (2000/60/EC) - Guidance Document No. 19: Guidance on Surface Water Chemical Monitoring Under The Water Framework Directive. Technical Report–2009- 025. URL. https://circabc.europa.eu/sd/a/e54e8583-faf5-478f-9b1141fda9e9c564/Guidance%20No%2019%20%20Surface%20water%20chemical%20m onitoring.pdf 27

615

Edmunds, W.M., Shand, P. 2008. Natural Groundwater Quality. Wiley-Blackwell , NewYork.

616 617 618

Guo, A.Z., and Bartram, J.K., Predictors of water quality in rural healthcare facilities in 14 low- and middle-income countries, J. Clean Prod. 237, 117836.

619

GWD, 2000. EU Groundwater Directive (2000/118/EC).

620 621 622

Gemitzi, A. 2012. Evaluating the anthropogenic impacts on groundwaters; a methodology based on the determination of natural background levels and threshold values. Environ. Earth Sci., 67, 2223-2237.

623 624 625 626

Hinsby, K., Condesso de Melo, M., Dahl, M. 2008. European case studies supporting the derivation of natural background levels and groundwater threshold values for the protection of dependent ecosystems and human health. Sci. Total Environ. 401(1-3), 1-20.

627 628 629

Hose, G.C. 2005. Assessing the Need for Groundwater Quality Guidelines for Pesticides Using the Species Sensitivity Distribution Approach. Hum. Ecol. Risk Assess. 11(5), 951–966.

630 631 632

Rodríguez-Huerta, E., Rosas-Casals, M., Hernández-Terrones, L.M., Water societal metabolism in the Yucatan Peninsula. The impact of climate change on the recharge of groundwater by 2030, J. Clean Prod. 235, 272-287.

633 634

Jenifer, M.A., Jha, M.K. 2018. Comprehensive risk assessment of groundwater contamination in aweathered hard-rock aquifer system of India. J. Clean Prod. 210, 853-868.

635 636 637 638 639 640

Karaaslan, Y. 2017. Development of a methodology for delineation of groundwater bodies: Gediz River Basin case study; In: Wolff, J. Führböter, J.F., Funk, P., Schleicher, K., Brinschwitz, D., Hartmann, A., Homrighausen, R., Lietzow, A., Meinert, M. (Eds), Zentralblatt für Geologie und Paläontologie, Teil I, Geologie; Special Issue: Groundwater protection between the conflicting priorities of sustainability and economy, VI, 232.

641 642 643

Mapoma, C., Xie, X., Pi, K., Liua, Y., Zhua, Y. 2016. Understanding arsenic mobilization using reactive transport modelling of groundwater hydrochemistry in the Datong basin study plot. Environ. Sci. Proces. Impacts 18(3), 371-385.

644 645

Matschullat, J., Ottenstein, R., Reimann, C. 2000. Geochemical background – can we calculate it?, Environ. Geol., 39(9), 990–1000.

646 647 648 649

MoFW (former Ministry of Forestry and Water Affairs of Turkey), 2017. Final Report by Fugro Sial Geosciences Consulting and Engineering Ltd. Co., Implementations for Determination and Assessment of Groundwater Quantity And Quality: Gediz Basin Pilot Study Project, Ankara, Turkey.

650 651 652

Molinari, A., Guadagnini, L., Marcaccio, M., Guadagnini, A. 2012. Natural back-ground levels and threshold values of chemical species in three large-scale groundwater bodies in Northern Italy. Sci. Total Environ. 425, 9-19. 28

653 654 655

Molinari, A., Guadagnini, L., Marcaccio, M., Guadagnini A., 2019. Geostatistical multimodel approach for the assessment of the spatial distribution of natural background concentrations in large-scale groundwater bodies, Water Res. 149, 522-532.

656 657 658 659 660

Müller, D., Blum, A., Hart, A., Hookey, J., Kunkel, R., Scheidleder, A., Tomlin, C., Wendland, F. 2006. D18: Final Proposal for a Methodology to Set up Groundwater Treshold Values in Europe. BRIDGE Project. Project co-funded by the European Commission within the Sixth Framework Programme (2002-2006), Contract no 006538 (SSPI).

661 662

Nakic, Z., Posavec, K., Bacani, A. 2007. A Visual Basic Spreadsheet Macro for Geochemical Background Analysis. Groundwater, 45(5), 642-647.

663 664 665 666

Nath, B.K., Chaliha, C., Bhuyan, B., Kalita, E., Baruah, D.C., Bhagabati, A.K. 2018. GIS mapping-based impact assessment of groundwater contamination by arsenic and other heavy metal contaminants in the Brahmaputra River valley: A water quality assessment study. J. Clean Prod. 201, 1001-1011.

667 668 669

Niu, J., Zhu, X., Parry, M. A. J., Kang, S., Ding, R., Environmental burdens of groundwater extraction for irrigation over an inland river basin in Northwest China, J. Clean Prod. 222, 182-192.

670 671 672

Papadopoulou Vrynioti , K. , Alexakis, D., Bathrellos, G.D. , Skilodimou, H.D., Vryniotis, D., Vassiliades, E. 2014. Environmental research and evaluation of agricultural soil of the Arta plain, western Hellas. J. Geochem. Explor., 136, 84-92.

673 674

Parrone, D., Ghergo, S., Preziosi, E. 2019. A multi-method approach for the assessment of natural background levels in groundwater. Sci.Total Environ. 659(1), 884-894.

675 676 677

Preziosi, E., Parrone, D., Del Bon, A., Ghergo, S. 2014. Natural background level assessment in groundwaters: probability plot versus pre-selection method. J. Geochem. Explor. 143, 43–53.

678 679

Ramírez, Y., Cisternas, L.A., Kraslawski, A., Application of House of Quality in assessment of seawater pretreatment technologies, J. Clean Prod. 148, 223-232.

680 681

Reimann, C., Filzmoser, P., Garrett, R. 2005. Background and threshold: critical comparison of methods of determination. Sci. Total Environ. 346(1-3),1-16.

682 683 684

Sellerino, M., Forte, G., Ducci, D. 2019. Identification of the natural background levels in the Phlaegrean fields groundwater body (Southern Italy), J. Geochem. Explor. 200, 181192.

685 686

SHW (State Hydraulic Works of Turkey), 2015. Hydrogeological Investigation Report for the Gediz Basin, Eser Project and Engineering, Ankara, Turkey.

687 688 689

Tedd, K., Coxon, K., Misstear, B., Daly, D., Craig, M., Mannix, A., and Williams, T.H. 2017. Assessing and Developing Natural Background Levels for Chemical Parameters in Irish Groundwater, Environmental Protection Agency, EPA Research Report, Ireland.

29

690 691

TBGW, 2012. The Turkish bylaw on the Protection of Groundwater against Pollution and Deterioration. Oficial Gazette Number: 28257.

692 693 694 695

Urresti-Estala, B., Carrasco-Cantos, F., Vadillo-Perez, I., Jimenez-Gavilan, P. 2013. Determination of background levels on water quality of groundwater bodies: A methodological proposal applied to a Mediterraean River Basin. J. Environ. Manag. 117, 121-130.

696 697 698

Vivona, R., Preziosi, E., Madé, B., Giuliano, G. 2007. Occurrence of minor toxic elements in volcanic-sedimentary aquifers: a case study in central Italy. Hydrogeol. J., 15, 11831196,

699 700 701

Wendland F, Hannappel S, Kunkel R, Schenk R, Voigt HJ, Wolter R. 2005. A procedure to define natural groundwater conditions of groundwater bodies in Germany. Water Sci Technol. 51(3–4), 249–57.

702 703 704 705

Wendland, F., Berhold, G., Blum, A., Elsass, P., Fritsche, J.-G., Kunkel, R., Wolter, R. 2008. Derivation of natural background levels and threshold values for groundwater bodies in the Upper Rhine Valley (France, Switzerland and Germany). Desalination 226(1-3), 160-168.

30