Abstract

Groundwater-monitoring network is a set of boreholes (wells) that is used to monitor the water table fluctuation and to detect groundwater contamination. In this research, the maximal covering location problem (MCLP) method is employed to discretize the area, and the LINGO modeling program is used to optimize the number of boreholes. Tabriz oil refinery at the northwest of Iran with high pollution potential that imposes a serious threat to the beneath multilayered aquifer was chosen to evaluate the feasibility of these techniques in field scale. The location and content of storage tanks, leakage history, groundwater flow direction, contaminated well location and the facilities leakage potential are considered as the weighting factors to calculate the number and location of the optimal wells. As a result of optimization, the initial estimated number of boreholes by the MCLP model for the study area is reduced from 349 to 184. A high density of optimal boreholes is allocated to refining zone and oil storing yard, especially near tanks containing dangerous substances due to their toxicity and potential for contaminating water. A vulnerability zoning map prepared using the analytical hierarchy process method indicates a suitable conformation between locations of the boreholes and the vulnerable areas.

HIGHLIGHTS

  • A multicriteria method is introduced to design and optimize the aquifer-monitoring network.

  • A reliable optimized monitoring network produces by the combination of MCLP and LINGO methods.

  • The initial number of 349 monitoring wells is optimized to 184 boreholes in an oil refinery.

  • A vulnerability zoning map prepared by the AHP method shows a suitable conformity with the density of boreholes.

Graphical Abstract

Graphical Abstract
Graphical Abstract

INTRODUCTION

Groundwater contamination by petroleum hydrocarbons is known as the most important problem that imposes threats to groundwater resources and humans (Al-Raoush et al. 2018). This problem is more serious in refineries or industrial petroleum areas, which causes the release of hydrocarbon contamination to underlying aquifers. Spilled or leaked petroleum products from storage tanks and pipes can accumulate in the subsurface as non-aqueous-phase liquids and lead to unacceptable risks for human health and ecosystems (Vaezihir et al. 2013). Therefore, groundwater as an essential and vital part of the water cycle needs careful protection and intelligent management against risks of petroleum hydrocarbon contamination released from refineries. Effective and efficient management of groundwater requires sufficient coverage of groundwater in both quantitative and qualitative aspects (Ayvaz & Elçi 2018). This goal can be achieved by establishing a groundwater-monitoring network that is designed to measure changes in groundwater quality or the groundwater level over time which plays an important role in aquifer restoration and prevention. Therefore, monitoring wells suspected to be polluted are of interest, and it is necessary to detect the pollutants as soon as possible.

It can be helpful to make informed management decisions regarding the status of complex groundwater systems by providing accurate and reliable data (Song et al. 2019). Groundwater-monitoring network is a set of wells that may consist of more or less randomly distributed wells drilled in the studied aquifer system (Ayvaz & Elçi 2018). Groundwater-monitoring networks are constructed to meet regulatory requirements, to monitor groundwater-level fluctuation and groundwater quality variations through time, to determine contamination sources and to the evaluation of the efficiency of remediation (Hassan 2003).

Data, which can be collected from a monitoring program, are the direction of groundwater flow, thickness of the unsaturated zone, monitoring of short- and long-term groundwater-level fluctuations, tracking of pollution caused by released contaminants and their fate and transport (Rosen 2003).

Typical monitoring networks consist of observation wells with a diameter of about 4–6 inches equipped with blind intubation (casing) in the unsaturated zone and with a screen (mesh) tube in the saturated zone. The design of a groundwater-monitoring network is defined as the selection of the monitoring borehole location and frequency (Ohmer et al. 2019). Since, the design and implementation of the monitoring network is very complex and expensive, optimal groundwater-monitoring networks is necessary to reduce monitoring effort, operation time and excessive costs. The number of monitoring wells depends on the objectives of the project and the methods of data analysis. A minimum number of monitoring wells with an optimum design that provides sufficient spatial coverage on groundwater quality becomes an important engineering optimization problem (Ayvaz & Elçi 2018), because increasing the number of monitoring wells more than optimized ones causes the extra cost and may not increase the accuracy of data significantly. Moreover, fewer monitoring wells will decrease the accuracy. Therefore, an ideal network is an optimized one from accuracy and financial points of view. The optimal design of a groundwater-monitoring network can be achieved by regarding several factors including hydrologic budget, geology of the aquifer, geography of the network area, potential sources, location of contaminants, the size of the monitoring network and incorporation of any regulatory requirements that may be needed to design (Rosen 2009).

Various statistical and modeling techniques are available to optimize and design groundwater-monitoring networks. The maximal covering location problem (MCLP), a numerical programming model, is one of the most common optimization model proposed in the literature that has been employed to find the best locations for public facilities in order to maximize the demand coverage (Church & ReVelle 1974).

The advantages of this method have been taken in various fields of studies such as redeployment of ambulance, locating first-aid facilities to help casualties in natural disasters and police patrol stations, location of light towers, cell phone towers and warning sirens and also design of emergency response systems.

For solving an MCLP, Linear, Integer, Nonlinear Programming and Global Optimization (LINGO), and algebraic modeling system, developed by the Linear, Integer, Nonlinear Programming and Discrete Optimizer (LINDO), can be used (Cunningham & Schrage 2004). Linear programming (LP) by making a few simplifying assumptions is known as one of the simplest ways to perform optimization that can be helpful to solve some very complex optimization problems (Vanderbei 2014).

For precise assessment of groundwater contamination, several methods have been investigated including the overlay and geographic information system (GIS)-based methods, the statistical methods and the process-based method (Mogaji 2018). The overlay and GIS-based methods are relatively simple and were commonly used in literature as decision support tools (Vaezihir & Tabarmayeh 2015; Hanssen et al. 2018; Machiwal et al. 2018; Richardson & Amankwatia 2018; Çelik 2019). The GIS-based analytical hierarchy process (AHP) is one of the most commonly used multicriteria evaluation methods developed by Saaty (1997) for decision-making in dealing with complex decision problems (Vu et al. 2019).

The design and optimization of groundwater-monitoring networks have widely been discussed in the literature over the last decades.

Vaezihir et al. (2012a, 2012b) presented a method to detect petroleum contaminants in the aquifer at the Shiraz oil refinery. Their work on source determination and modeling of fate and transport of BTEX (benzene, toluene, ethylbenzene and xylenes) contaminants revealed that the current monitoring network needed to be augmented. As a result, two strategies were defined for site clean-up: (1) removal of BTEX from the entire site down to MCL and (2) remediation of the site to a level that controls the migration of the BTEX compounds before they reach the nearest discharge wells. Farlin et al. (2019), by measuring contaminant concentration and taking solute transit time into account, employed a method to select and optimize groundwater quality-monitoring network stations from a larger choice of fixed locations such as springs or existing observation wells. In the presented method, at first, the empirical distribution of selected solutes (contaminants and age tracers) was calculated from a large dataset consisting of as many sampling points as possible and compared to the distributions obtained solely from the stations of the existing or planned monitoring network. In the case of high differences, the network was modified by replacing some of stations or adding new ones. At the second step, the distribution of the optimized network was compared to the total distribution of other contaminants and groundwater age tracer. The use of the presented method was applied for the main aquifer of the country of Luxembourg, where it is indicated that the representativeness of the available monitoring network can be improved by replacing two stations.

Hudak & Loaiciga (1993) presented an approach for designing a detection-based groundwater quality-monitoring network in multilayered groundwater flow systems. The presented methodology in this study is a viable approach to design a detection-based groundwater quality-monitoring network in multilayered groundwater flow systems. Susceptibility to contamination at points throughout a model domain was quantified by weight values. The weights were utilized in a mathematical programming model, which was selected monitoring sites in each of several hydrostratigraphic intervals comprising an overall model domain. The approach was tailored to the early detection of a contaminant release as opposed to the characterization of a contaminant plume. Maqsood et al. (2004) designed programs for groundwater monitoring in an area contaminated by oil in western Canada based on the behavior of contaminant transport in groundwater and the optimization of a monitoring network using statistical trend analysis. Bierkens (2006) specified the number of monitoring wells with the minimum total cost in an industrial area using a combination of stochastic simulation and a cost model. Thakur (2015) identified and investigated new approaches by using statistical and geostatistical methods for groundwater-monitoring network optimization to improve groundwater-monitoring strategies. The study has demonstrated that the existing monitoring network could agreeably be optimized using the presented statistical and geostatistical methods without the obvious fear of losing essential information from the current monitoring network. Hosseini & Kerachian (2017) proposed a new methodology to improve the reliability of groundwater-monitoring networks using combined numerical, geostatistical and neural network-based simulation models. The results showed that a new sampling configuration with 35 and 7 monitoring stations leads to a more efficient monitoring network than the existing one containing 52 monitoring stations. The mean variance estimation errors of all scenarios decrease significantly when compared with that of the existing monitoring network. Mirzaie-Nodoushan et al. (2017) proposed a method which relied on time series of groundwater levels to optimal design of groundwater-level monitoring networks in the Eshtehard aquifer that showed successful results in achieving accuracy and cost-effectiveness. This study relied on time series of groundwater levels to design groundwater-monitoring networks. The optimization algorithm employed in this paper considered as the entire area of the aquifer in search of the best monitoring sites. Ayvaz & Elçi (2018) presented the identification of the optimum groundwater quality-monitoring network using the genetic algorithm (GA)-based optimization approach. Findings indicated that the proposed approach significantly reduces the number of monitoring wells with a relatively small deviation of the spatial distribution of the WQI values.

Badri & Sedaghat (2017) applied the GA and AHP to provide oil spill risk map in the Persian Gulf. GA was used to obtain the optimum placement of the vessels for several oil spill scenarios. Meanwhile, by identifying and analyzing six major factors affecting the damage of ecosystem, the risk map of oil spill was generated for the Persian Gulf.

The objective of this research is the application of a MCLP and LINGO modeling program to design and optimize groundwater-monitoring network for detecting groundwater contamination and monitoring of groundwater-level fluctuation in Tabriz Oil Refinery and Petrochemical Complex (TRPC). Tabriz oil refinery began operations in 1977. This area was attacked by Iraqi aircrafts 21 times throughout the Iran–Iraq war (1980–1988). Due to these attacks and corrosion of the facilities, it is very possible that oil materials may leak into the aquifer beneath. In addition, Tabriz oil refinery is producing 80,000 bpd barrels per day of oil products such as regular and ATK, fuel oil, diesel oil and premium gasoline. In this industrial area, a variety of oil products and related hydrocarbon compounds can be released into groundwater due to various potential sources of leakage. If leakage occurs from processing units, storage tanks and conveyance lines of oil products, it is possible that dangerous oil pollution penetrates into the groundwater. Therefore, groundwater resources face the increasing risk of contamination. As a result, due to the lack of groundwater-monitoring network, there is no available reliable and comprehensive information about possible contamination. So, design and installation of a monitoring network in this industrial area is a necessity. In this research for the first time, a multicriterion method using the combination of MCLP and LINGO modeling program is applied to design and optimize groundwater-monitoring network in the oil refinery site, and the groundwater contamination potential map was provided using analytical methods to evaluate the accuracy of the designed network.

MATERIALS AND METHODS

Study area

The TRPC is located in the northwest of Iran, more specifically in the East Azerbaijan province, and adjacent to the Tabriz city (Figure 1). The Tabriz oil refinery was founded in 1977 with a refining capacity of 115,000 barrels per day of crude oil, which, by implementing development plans, has been enlarged to 80,000 bpd at the present time. The feedstock is crude oil and the main products include gasoil, gasoline, fuel oil, kerosene, naphtha, sulfur, ethane, vacuum bottom and liquid petroleum gas (LPG). The Tabriz petrochemical company (TPC) is located next to the refinery in about 380 ha. TPC is a producer of raw polymers such as polyethylene, polystyrene and ABS, each in different grades. Main feed and consumed raw materials are naptha, LPG, acrylonitrile, mineral oil, polybutadiene (PBR) and other chemicals.

Figure 1

Location of the study area.

Figure 1

Location of the study area.

The alluvial aquifer in the study area consists of at least two layers, of which the shallow layer is salty and not suitable for use. The second layer is confined and has fresh water. Because any leakage from the industries that might threaten the lower aquifer needs first to pass through the upper aquifer, in this research, the monitoring network is designed for the upper, unconfined aquifer. However, flow between the two aquifers is possible due to wells that penetrate both; thus, it is possible that industrial, agricultural and domestic contaminants that reach the upper aquifer could contaminate the confined aquifer beneath it as well (Yekom Consulting Engineers 2012). The general direction of groundwater flow in the area is from the southeast to the west and northwest (Figure 2).

Figure 2

Pumping wells of the study area (circled dots show the sampling wells).

Figure 2

Pumping wells of the study area (circled dots show the sampling wells).

Groundwater sampling and analysis

The monitoring network should be denser in polluted areas, so it is necessary to investigate any pollution of the aquifer based on sampling from existing wells and boreholes. Therefore, 13 groundwater samples were collected in May 2014 and were analyzed for TPH (total petroleum hydrocarbon) using the GC–MS (gas chromatography with mass spectrometry detector) method (Figure 2). Samples were collected in 60 mL amber glass vials and 0.6 mL of sodium azide (NaN3) was added to each sample to prevent biodegradation of hydrocarbons.

The standard static headspace method was used for the extraction of the analytes. The headspace gas was injected directly into the inlet of GC–MS (GC Agilent 6890N, MS Agilent 5973N, USA) with a 30 m, 0.250 mm I.D., 0.25 μm coated column. The carrier gas was helium with a flow rate of 2 mL/min. The chromatographic conditions were injection port temperature −275 °C, initial column temperature −40 °C for 5 min, heating rate −7 °C/min to a final temperature of 180 °C after 10 min and MS detector temperature −325 °C. The results of analysis by GC–MS were obtained as fingerprint diagrams for each sample.

Groundwater-monitoring network

Meyer & Brill (1988) presented a method for the MCLP to insert a well in a monitoring network under conditions of uncertainty and heterogeneity. This is a programming formulation for the design and optimization of the location of a fixed number of points in the network of demanded nodes. For example, in the case of the addition of a number of firefighting centers in a city to an existing network, the city map is divided into a set of demanded nodes, including a large number of possible sites to put new facilities for each node has a weight which indicates related population in that location, and the firefighting centers allocated to the closest node. For emergencies, it is not reasonable for a demand node to be located more than a threshold distance from a firefighting center; in this case, it is considered an ‘uncovered node’. This threshold value is called the ‘maximum service distance’. If at least one center exists within its maximum service distance, it is considered as a covered node. If facilities are insufficient for reaching the entire coverage at the maximum service distance, it is necessary to look for more nodes to cover these facilities. In the present research, this theory was employed to determine the optimum number of monitoring boreholes to cover the entire site without posing unacceptable errors in the accuracy of the monitoring program.

To optimize the borehole number in the primary network, first an objective function should be defined for the limitation of the constraints. In this function, a weight is assigned to each node (wi) and each node is covered (yi) by the surrounding nodes, which is defined as the maximum service distance for the nodes. Where the node is covered by the surrounding ones, it is assigned a value of 1; otherwise, it is considered equal to zero. All nodes should be checked for coverage mathematical formulation of the objective function. Equation (1) and constraints of Equation (2) for this problem are as follows (Church & ReVelle 1974):
formula
(1)
formula
(2)
formula
(3)
formula
(4)
formula
(5)
formula
(6)
where
  • wi is the population at node i, I is the set of demand nodes in the discretized network, J is the set of prospective nodes for siting facilities, Jp is the set of nodes j occupied by pre-existing facilities.

  • Ni:;

  • dij is the shortest distance from node i to node j; S is the cover distance threshold (maximal service distance); xj equals to 1 if a facility is located at site j; 0 otherwise; yi equals to 1 if node i is covered; 0 otherwise; P is the total number of facilities to be located.

LINGO is a mathematical modeling language and optimizer program that supports the analysis and solution of a broad range of mathematical programming. LINGO is widely used to solve optimization problems including linear, nonlinear, integer and quadratic programming (Sharif & Swamy 2014).

One of the simplest ways to perform optimization and the earliest formulated problems in mathematical programming is LP which is dealing with the optimization of a linear objective function under linear equality and linear inequality (Vanderbei 2014). Index, variable, constraint and objective function are the components that must be considered in formulating LP (McCarl & Spreen 1997). Algebraically, LP can be formulated in the following form (Vanderbei 2014):
formula
(7)
formula
(8)
formula
(9)
where x1 …, xn are called variables of the LP; n and m are the number of variables and constraints, not including the negativity constraints, respectively.

In this research, discretization of the study area was carried out using the MCPL method and then the boreholes were optimized through the LINGO modeling system by application of the LP technique.

Assigning weights to the nodes

One of the main aims of optimization is to find the optimum frequency and locations of the monitoring nodes by finding factors that achieve the highest performance of monitoring wells. Desirable designing and optimization of monitoring network is attainable by considering the affected multiple criteria on groundwater quality and quantity and also by determining the weights of the criteria through considering its importance in pollution by the possible source to each node.

In this research, the location and contents of aboveground storage tanks (ASTs), leakage potential of the facilities, leakage history, contaminated points and groundwater flow direction importance were identified and weighted as the main affecting groundwater criteria.

Among the hydrocarbon components, BTEX and MTBE (methyl tert-butyl ether) with a high molar mass of aromatic components are the most important detected petroleum products in groundwater that are highly soluble and mobile in water and have potential threats to human health (EPA 2015; Stefanakis 2020). MTBE is an additive and is used as an oxygenating factor to enhance the gasoline octane number (AnishRaman et al. 2014). MTBE is a synthetic chemical compound that is derived from the catalytic reaction of methanol and isobutylene (Axelrod & Coleman 2019). BTEX is generally found in crude oil and its refined products are especially found in gasoline with an average weight of 18% and 26.93 molar percent (Table 1). BTEX is relatively soluble in water but, because of their absorption by sediments, the speed of BTEX movement is slower than MTBE in groundwater. MTBE is more soluble than another gasoline compounds like benzene, has less tendency to be absorbed by sediments and is more resistant to biodegradation than BTEX (EPA 1998).

Table 1

Molar percent of BTEX in typical oil materials produced in Iran (adapted from Vaezihir et al. 2012a, 2012b)

Type of oil materialsBTEX molar percent
Gasoline 26.93 
ATK 6.01 
Kerosene 4.09 
Crude oil 3.24 
MTBE 0.72 
Gas oil 0.63 
Fuel oil 0.19 
Type of oil materialsBTEX molar percent
Gasoline 26.93 
ATK 6.01 
Kerosene 4.09 
Crude oil 3.24 
MTBE 0.72 
Gas oil 0.63 
Fuel oil 0.19 

Therefore, it is important to consider these factors in designing groundwater-monitoring networks, which show the importance of tanks and other facilities in the case of leakage.

The analytical hierarchy process

The analytic hierarchy process (AHP) developed by Saaty (1977) is a theory for dealing with complex technological, economical, and sociopolitical problems. This is an attempt to simulate the modeling of problems, in which each problem has its own specific model and terminology (Saaty & Alexander 1989). The AHP method is used to derive a scale of relative importance of alternatives or criteria from pairwise comparisons (Richardson & Amankwatia 2018). This method can be useful for evaluating and weighting criteria in multicriterion decision-making to reach an acceptable solution. This method is based on the pairwise comparison matrix and requires the active participation of decision-makers in reaching an agreement decision (Dey 2004).

Pairwise comparison matrix

At the last step of AHP, data of the criteria are presented to a number of experts who will carry out a pairwise comparison of the criteria with respect to each risk factor. The outcome of the comparison is a matrix that ranks the criteria in order of likelihood of contamination. Experts need to rank each factor using a scale of 1–9. For example, if there are two criteria for a level of risk, then pairwise comparison score is regarded as 1. A score 9 will be allocated if there is a much stronger criterion than the other (Dawotola et al. 2010). After creation of binary comparison matrices, relative weight is derived for different criteria. The relative weights of the criteria for each level with respect to a higher-level element are calculated as the components of the element associated with the largest eigenvalue of their comparison matrix. Then, the weights will be combined through the hierarchy method to produce composite weights (Dey 2004).

Weighted overlay method

In many suitable analyses, some eligible criteria are more important than others are. Often, it is the expectation in a location search to be able to compare several suitable candidates whether, and how strongly, they meet a set of criteria by differing importance. Using the layer principle, one can easily extend the overlay by assigning levels of importance to each criterion. A numerical weighting factor is assigned to each thematic layer according to its relative importance compared to all other layers. After that, the weighted layers are covered with the previous layers. This is called weighted overlay (Flitter et al. 2013). In this research, weighted overlay is performed for selected criteria including contamination points (sampling points), areas with a history of leakage, land use, soil permeability and flow direction that were rated and weighted based on the AHP method results.

Linear combination of all factors can be written as the following equation:
formula
(10)
where S is the total index, wi is the weight of factor i and xi is the rating score of factor i. According to the results of this method, the area can be divided into classes based on the research objective.

RESULTS AND DISCUSSION

Investigation of groundwater pollution by petroleum products

According to the results of hydrocarbon component analyses, 2 out of 13 samples indicated significant contamination (Figure 3 and Table 2). Chromatograph fingerprints of contaminated samples (W2 and W7) are presented in Figure 3(a) and 3(b), which show a high abundance of hydrocarbon compounds, compared to a background sample in Figure 3(c). It implied that this area of the groundwater might be contaminated by leakage or spill of oil materials. In the monitoring network design process, the results of this sampling were used to assign a higher weight to the contaminated areas. It meant that a dense network was needed in the contaminated part of groundwater.

Table 2

The results of TPH analysis obtained by chromatograph (Figure 3)

SamplePeak no.Registration time (min)Hydrocarbon materialPeak heightArea below the peakMaximum of peak height (%)The amount of hydrocarbon material, volume percentage (%)
W2 17.61 Phenol, 5-methyl-2-(1-methylethyl) 74,513 15,710,544 100 37.38 
18.21 Phenol, 5-methyl-2-(1-methylethyl) 29,594 812,108 5.17 1.93 
18.29 1-(P-methoxybenzamido)-2-(cyclohexene-1-yl) 33,425 1,025,771 6.53 2.44 
18.49 Cyclohexasiloxane, dodecamethyl 118,974 12,132,788 77.23 28.87 
18.79 2- Cyclohexen-1-one, phenol, 5-methyl- 2-(1-methylethyl) 29,578 928,945 5.91 2.21 
20.15 2,6-Octadien-1-ol, 3,7-dimethyl 23,428 1,524,736 9.71 3.63 
11 22.09 Cycloheptasiloxane, tetradecamethyl 169,708 5,643,436 35.92 13.43 
16 25.02 [[4-[1,2-Bis[(trimethylsilyl)oxy]ethyl]-1,2-phenylene]bis(oxy)]bis(trimethylsilane) 58,145 2,360,736 0.47 0.38 
W7 14 20.31 Tetradecene 40,545 5,551,328 82.44 18.11 
30 23.47 1-Pentadecene, 2-methyl-56 23,075 1,042,594 15.48 3.4 
31 23.8 2-Tetradecene, 1-Hexadecene 139,356 6,733,545 100 21.97 
32 23.9 7-Hexadecene 46,539 3,631,995 53.94 11.85 
SamplePeak no.Registration time (min)Hydrocarbon materialPeak heightArea below the peakMaximum of peak height (%)The amount of hydrocarbon material, volume percentage (%)
W2 17.61 Phenol, 5-methyl-2-(1-methylethyl) 74,513 15,710,544 100 37.38 
18.21 Phenol, 5-methyl-2-(1-methylethyl) 29,594 812,108 5.17 1.93 
18.29 1-(P-methoxybenzamido)-2-(cyclohexene-1-yl) 33,425 1,025,771 6.53 2.44 
18.49 Cyclohexasiloxane, dodecamethyl 118,974 12,132,788 77.23 28.87 
18.79 2- Cyclohexen-1-one, phenol, 5-methyl- 2-(1-methylethyl) 29,578 928,945 5.91 2.21 
20.15 2,6-Octadien-1-ol, 3,7-dimethyl 23,428 1,524,736 9.71 3.63 
11 22.09 Cycloheptasiloxane, tetradecamethyl 169,708 5,643,436 35.92 13.43 
16 25.02 [[4-[1,2-Bis[(trimethylsilyl)oxy]ethyl]-1,2-phenylene]bis(oxy)]bis(trimethylsilane) 58,145 2,360,736 0.47 0.38 
W7 14 20.31 Tetradecene 40,545 5,551,328 82.44 18.11 
30 23.47 1-Pentadecene, 2-methyl-56 23,075 1,042,594 15.48 3.4 
31 23.8 2-Tetradecene, 1-Hexadecene 139,356 6,733,545 100 21.97 
32 23.9 7-Hexadecene 46,539 3,631,995 53.94 11.85 
Figure 3

Chromatograph (fingerprint) of TPH for groundwater samples: (a) W2, (b) W7 and (c) W13 (as background).

Figure 3

Chromatograph (fingerprint) of TPH for groundwater samples: (a) W2, (b) W7 and (c) W13 (as background).

The highest retention times of total hydrocarbon materials in the sample W2 are phenol, 5-methyl-2-(1-methylethyl), cyclohexasiloxane, dodecamethyl and cycloheptasiloxane, and tetradecamethyl compounds with 37.38, 28.87 and 13.43%, respectively. In decreasing order, 2-tetradecene, 1-hexadecene, tetradecene and 7-hexadecene are the most common hydrocarbon materials in sample W7 with 21.97, 18.11 and 11.85%, respectively.

Designing the monitoring network

Designing the groundwater-monitoring network was carried out in following steps.

Discretization of the study area

The study area was discretized into 349 nodes using the MCLP method, in which the location of nodes represents the sites for boreholes. Cell size (and thus, node density) was assigned based on the importance of the location, which was defined by risk of contamination in a given part of the refinery. ASTs, the refining area, transport pipelines, loading areas and areas with a history of leakage were regarded as the important areas for high risk of contamination. For example, the location of tanks was discretized with a nodal separation of 100 m, whereas building areas were discretized with lower density (400 m; Figure 4). The number of these nodes as the location of hypothetical boreholes needs to be optimized using the LP method.

Figure 4

Map of the distribution of nodes in the study area.

Figure 4

Map of the distribution of nodes in the study area.

Optimization of the groundwater-monitoring network

To optimize the primary assigned nodes (discretized area), it is necessary to determine the objective function and constraints. The objective function is optimizing the model to the number of nodes with the maximum possible accuracy by considering its distance to refinery facilities like tanks or office site and their relative position to the groundwater flow direction and minimum cost that is estimated by calculating the required budget for the monitoring process such as wells drilling and maintenance, sample analysis and so on. The general objective function is defined as follows:
formula
(11)
where yi, wi and xi are the covered node, the weight assigned to each node and the initial node location, respectively.

The following steps were carried out to develop groundwater-monitoring network by determining the objective function under constrains.

Numerical code development

A matrix was prepared to determine the number of nodes, and then the maximum service distance and coordinates of the nodes were defined for different zones. The minimum distance between each node relative to surrounding nodes was determined. If the distance between a node and at least one of the neighbor nodes is less than or equal to the maximum service distance, it means that this node has been covered by the surrounding nodes. These constraints were determined by writing a separate MATLAB program for each of the zones.

Determination of the constraints

The numerical code needs to evaluate the number of nodes and, at the same time, it should satisfy the constraints. The main constraint is that any point in the study area must be covered by at least one node, requiring that all places in the TPRC area should be covered by one monitoring well established by the MCLP method, and that any release of contaminants into groundwater can be detected by at least one of the surrounding monitoring wells.

Selection and weighting of criteria

The most important step in this method is to assign the weight to each node. In this study, contaminated points, groundwater flow direction, storage tank location, storage tank content, leakage potential of the facilities and leakage history were selected as the main factors for weighting each node. Also, different weights were assigned to the nodes according to the importance of the petroleum materials based on their potential risk to the environment (Table 3). Due to the high risk of mono-aromatic compounds, the BTEX content of the oil products was adopted as the scale for weight assignment.

Table 3

Weights assigned to the area based on the tank contents

Tank contentWeight
Benzene, toluene, ethylbenzene 0.26 
Gasoline 0.21 
Crude oil, kerosene, naphta 0.16 
Solvent, condensate 0.13 
Gas oil 0.11 
ATK 0.08 
Fuel oil, pitch, isophide 0.05 
Tank contentWeight
Benzene, toluene, ethylbenzene 0.26 
Gasoline 0.21 
Crude oil, kerosene, naphta 0.16 
Solvent, condensate 0.13 
Gas oil 0.11 
ATK 0.08 
Fuel oil, pitch, isophide 0.05 

Running optimization model

The weights and derived constraints were used to solve the optimization objective with the LINGO program in order to obtain the number of boreholes required in each region. The total number of proposed optimal boreholes was 184, of which 29, 73 and 57 were assigned to the distribution unit, refinery unit and petrochemical unit, respectively. The remaining 25 boreholes were needed to fill in the areas surrounding these units (Figure 5).

Figure 5

Final design for the groundwater-monitoring network.

Figure 5

Final design for the groundwater-monitoring network.

Contamination potential map of groundwater

The AHP model along with the weighted overlay method was employed to prepare zonation map that indicates the areas with high probability of contamination. In order to prioritize influencing factors on groundwater contamination, the pairwise comparison matrix was extracted using expert assessment for five parameters including contaminated points, areas with a history of leakage and bombardment, land use (such as loading areas, refining units and reservoirs), soil permeability and flow direction (Table 4). The final weights and rates were computed for each criterion based on the result of pairwise comparison and the consistency ratio between the expert's opinions (Table 5).

Table 4

Pairwise comparison matrix for multiple criteria

CriterionContamination of sampling pointsLeakage and bombarded areasLand useGroundwater flow directionWeight of relative
Contamination of sampling points 0.47 
Leakage and bombarded areas 0.5 0.27 
Land use 0.33 0.5 0.16 
Groundwater flow direction 0.25 0.33 0.5 0.10 
CriterionContamination of sampling pointsLeakage and bombarded areasLand useGroundwater flow directionWeight of relative
Contamination of sampling points 0.47 
Leakage and bombarded areas 0.5 0.27 
Land use 0.33 0.5 0.16 
Groundwater flow direction 0.25 0.33 0.5 0.10 
Table 5

The rates and weights assigned to each parameter based on their relative importance on groundwater contamination

ParameterRateWeight
Contaminated points  
 High contamination 
 Moderate contamination 
 Low contamination 
 Very low contamination 
 Not detected contamination 
Land use  
 Reservoirs and refining units 
 Facilities 
 Other areas 
Leakage history and bombardment impact  
 Underground installations, dumping of hydrocarbon wastes 
 No leakage 
Groundwater flow direction (m)  
  <1,310 
 1,310–1,315 
 1,315–1,333 
 1,333–1,342 
 1,338–1,342 
 1,342 < 
ParameterRateWeight
Contaminated points  
 High contamination 
 Moderate contamination 
 Low contamination 
 Very low contamination 
 Not detected contamination 
Land use  
 Reservoirs and refining units 
 Facilities 
 Other areas 
Leakage history and bombardment impact  
 Underground installations, dumping of hydrocarbon wastes 
 No leakage 
Groundwater flow direction (m)  
  <1,310 
 1,310–1,315 
 1,315–1,333 
 1,333–1,342 
 1,338–1,342 
 1,342 < 

The groundwater contamination potential map was then generated following the weighted linear combination method. The preparation of the criterion maps is described as below.

Contaminated points

Results of sample analysis are considered as the most important parameter to produce contamination risk maps. According to sampling and analysis results, the points located inside and downstream of the refinery and two points inside the petrochemical complex showed more contamination. Therefore, this area as an effective factor in groundwater contamination allocated high rates relative to the other points (Figure 6(a) and Tables 4 and 5).

Figure 6

Weighted map of groundwater potential for contamination based on (a) TPH, (b) leak history, (c) land use and (d) groundwater flow direction.

Figure 6

Weighted map of groundwater potential for contamination based on (a) TPH, (b) leak history, (c) land use and (d) groundwater flow direction.

Land use

Land use in an oil refinery means the type of facilities and installation that may release pollution to the environment. Accordingly, the loading site, location of the reservoirs (tank yard) and refining units are considered as areas with high potential of oil release (Figure 6(b) and Tables 4 and 5).

Leakage history and bombardment impact

Pollution may take place as leakage from underground installations or by the dumping of hydrocarbon wastes or the accidental release of refined product. As mentioned above, Tabriz Refinery has been subjected to bombardment by Iraqi air forces 21 times. So, the release of oil contaminant to the environment has the main role in the pollution of soil and groundwater of the site. On the other hand, due to age of the installation of the refinery (opened at 1978), the underlying conditions that govern the progress of corrosion must be considered as potential sources of groundwater contamination. Spilled or leaked light hydrocarbons will move downwards through soil until they reach the water table where they will float as a lens over saturated groundwater. Therefore, points with leakage and bombarded areas are considered as areas with high vulnerability (Figure 6(c) and Tables 4 and 5).

Groundwater flow direction

Depending on its physical, chemical and biological properties, a contaminant that has been released into the environment may move through an aquifer under government of advection as the main phenomenon of mass transport. Because of this slow movement, contaminants tend to remain concentrated in the form of a plume (EPA 2015).

Theoretical forecasts indicate that contaminated groundwater may extend in the flow direction from less than one meter to several thousand meters from the site of leakage, depending on the rate of sorption, biodegradation and the water velocity through the porous media (Duffy et al. 2003).

Groundwater flow at the study area occurs from the southeast to the west, and so, the potential of a location for contamination is considered from less to more in the flow direction (Figure 6(d) and Tables 4 and 5).

Soil permeability

Groundwater is subjected to contamination unless protected by a low permeability layer such as clay. Organic contaminants such as halogenated organic compounds, hydrocarbons or other organic compounds, which contaminate the soil matrix, can be resistant for a long time (Abdel-Moghny et al. 2012). Soils that are porous and permeable tend to transmit water and certain types of contaminants easily to the aquifer (EPA 2015). Therefore, soil permeability is an important factor influencing groundwater contamination. The soil permeability of Tabriz Refinery area is evaluated by its lithology based on borehole log data, which revealed that it consists of silt and clay. Due to the presence of silt and clay in surface layers, low vulnerability can be assigned to the site in terms of the infiltration of ground surface-leaked contaminants into groundwater. Due to the uniformity of the soil texture in the whole region, this criterion was not considered in the zonation process.

A groundwater contamination map was obtained by combining classified data layers in the GIS environment using the overlaying technique based on designed weight to each criterion in Table 5. As shown in Figure 7, the resultant Tabriz oil refinery groundwater contamination map was categorized into various regions from low to high vulnerability, indicating high vulnerability of the region in the central and southwest parts of the refinery (Figure 7).

Figure 7

Contamination risk map and compatibility with the designed monitoring network.

Figure 7

Contamination risk map and compatibility with the designed monitoring network.

In order to adjust the monitoring network with the results of contamination vulnerability zoning map, the result of the MCLP model was projected on the AHP risk map. The map shows a suitable conformity between the location of the boreholes and the vulnerable areas (Figure 7). Based on this assessment, about 75% of the proposed boreholes has been located in the vulnerable parts of the study area, indicating good computability between two methods.

CONCLUSIONS

This study introduced a multicriteria analysis to design and optimize the groundwater-monitoring network using MCLP and LINGO modeling programs in the case of a petroleum industrial area with various potential sources of pollution. Tabriz oil refinery and petrochemical and distribution units located on top of the multilayered aquifer imposes a high risk of leakage into the groundwater and is a serious threat to groundwater quality. Based on the groundwater hydrocarbon compound analysis, 2 out of 13 samples that belong to the boreholes located inside the refinery and downstream of it displayed high concentrations of oil products in the groundwater. Therefore, in order to restore and prevent groundwater from pollution, the MCLP model and LINGO modeling program were employed for design and optimization of initial borehole numbers in monitoring network. The study area was discretized into cells and nodes, and the constraints of the optimization were determined. Then, weights were assigned to the nodes according to the risk of pollution at their locations. The optimal number of boreholes was estimated to solve the objective function and satisfy constraints. The initial number of the boreholes was 349, which was reduced to 184 by optimization. According to the results, optimization decreased the number of boreholes up to 52.72%, with greatest density in the refining zone and oil storing, especially tanks containing dangerous substances due to their toxicity and potential for water contamination. A lower density of monitoring wells was proposed for areas outside of the refinery and petrochemical boundaries. A vulnerability zoning map, which was prepared using the AHP method and by combining five considered parameters including contaminated points, areas with a history of leakage and bombardment, land use (such as loading areas, refining units and reservoirs) and soil permeability, indicated a suitable conformity between the location of the boreholes and the vulnerable areas, which indicates the fair optimization and designing of the monitoring network.

The method introduced in this research is a low cost optimization procedure that the required parameters can be obtained easily. In addition, in this method, the basis on which the criteria are selected and weighted is generally simple and understandable. Therefore, these methods can be applicable and economical methods in industrial sites, especially at industrial sites with various potential sources of contamination.

ACKNOWLEDGEMENTS

This work was supported by Tabriz Oil Refining Company and University of Tabriz. The authors express their gratitude for the careful reviews to Derek C. Ford from McMaster University because he raised many important points, which improved the original manuscript.

DATA AVAILABILITY STATEMENT

Some of the data, models, or codes that support the findings of this study are available from the corresponding author upon reasonable request.

REFERENCES

Abdel-Moghny
T.
Mohamed
R. S.
El-Sayed
E.
Mohammed Aly
S.
Snousy
M. G.
2012
Effect of soil texture on remediation of hydrocarbons-contaminated soil at El-Minia district, Upper Egypt.
ISRN Chemical Engineering
2012
,
406598
.
Al-Raoush
R.
Ngueleu
S.
Rezanezhad
F.
Van Cappellen
P.
2018
Groundwater pollution by petroleum derived contaminants in coastal semiarid environment
. In
Qatar, Doha
,
Foundation Annual Research Conference Proceedings (Vol. 2018, No. 1, p. EEPD709). HBKU Press, Hamad bin Khalifa University Press
.
AnishRaman
C.
Varatharajan
K.
Abinesh
P.
Venkatachalapathi
N.
2014
Analysis of MTBE as an oxygenate additive to gasoline
.
International Journal of Engineering Research and Applications
4
,
712
718
.
Axelrod
M. G.
Coleman
S. T.
2019
U.S. Patent No. 10,364,204
.
U.S. Patent and Trademark Office
,
Washington, DC
.
Badri
M. A.
Sedaghat
A.
2017
Application of Genetic Algorithm and Analytic Hierarchy Process to Generate an Oil Spill Risk Map
.
Research Institute for Subsea Science & Technology, Isfahan University of Technology
,
Iran
.
ISSN 2378-3184
.
Bierkens
M. F. P.
2006
Designing a monitoring network for detecting groundwater pollution with stochastic simulation and a cost model
.
Stochastic Environmental Research and Risk Assessment
20
(
5
),
335
351
.
Church
R.
ReVelle
C.
1974
The maximal covering location problem
.
Papers of the Regional Science Association
,
32
(
1
),
101
118
.
Cunningham
K.
Schrage
L.
2004
The LINGO algebraic modeling language
. In
Modeling Languages in Mathematical Optimization
(
Kallrath
J.
, ed.).
Springer
,
Boston, MA
, pp.
159
171
.
Dawotola
A. W.
Van Gelder
P. H. A. J. M.
Vrijling
J. K.
2010
Multi Criteria Decision Analysis framework for risk management of oil and gas pipelines.
In
Reliability, Risk and Safety
(
Ale
B. J. M.
Papazoglou
I. A.
Zio
E.
, eds).
CRC Press
,
London
, pp.
307
314
.
Dey
P. K.
2004
Decision support system for inspection and maintenance: a case study of oil pipelines
.
IEEE Transactions on Engineering Management
51
(
1
),
47
56
.
Duffy
J. J.
Peake
E.
Mohtadi
M. F.
2003
Oil spills on land as potential sources of groundwater contamination
.
Environment International Journal
3
,
107
120
.
EPA
1998
Remediation of MTBE Contaminated Soil and Groundwater.
Environmental Protection Agency, United State, EPA. 510-F-97-015
.
EPA
2015
Ground Water Contamination, Wellhead Protection: A Guide for Small Communities.
Office of Research and Development Office of Water, Washington, DC. Chapter 3. EPA/625/R-93/002
.
Flitter
H.
Laube
P.
Lüscher
P.
Rogers
S.
Hägi
S.
2013
Suitability analysis.
Geographic Information Technology Training Alliance (GITTA)
,
1
25
.
Hanssen
F.
May
R.
van Dijk
J.
Rød
J. K.
2018
Spatial multi-criteria decision analysis tool suite for consensus-based siting of renewable energy structures
.
Journal of Environmental Assessment Policy and Management
20
(
03
),
1840003
.
Hassan
A. E.
2003
Long-Term Monitoring Plan for the Central Nevada Test Area (No. DOE/NV/13609-30; DRI Pub. No. 45201)
.
Nevada Site Office
,
Las Vegas, NV
,
USA
.
Machiwal
D.
Jha
M. K.
Singh
V. P.
Mohan
C.
2018
Assessment and mapping of groundwater vulnerability to pollution: current status and challenges
.
Earth-Science Reviews
185
,
901
927
.
McCarl
B. A.
Spreen
T. H.
1997
Applied Mathematical Programming Using Algebraic Systems
.
Texas A&M University
,
Cambridge, MA
.
Mirzaie-Nodoushan
F.
Bozorg-Haddad
O.
Loáiciga
H. A.
2017
Optimal design of groundwater-level monitoring networks
.
Journal of Hydroinformatics
19
(
6
),
920
929
.
Ohmer
M.
Liesch
T.
Goldscheider
N.
2019
On the optimal spatial design for groundwater level monitoring networks
.
Water Resources Research
55
(
11
),
9454
9473
.
Richardson
C. P.
Amankwatia
K.
2018
GIS-based analytic hierarchy process approach to watershed vulnerability in Bernalillo County, New Mexico
.
Journal of Hydrologic Engineering
23
(
5
),
04018010
.
Rosen
M. R.
2003
Groundwater monitoring networks
. In:
Groundwater
(
Silveira
L.
ed.).
Encyclopedia of Life Support Systems (EOLSS), Eolss Publishers
,
Oxford
,
UK
, pp.
2
.
Rosen
M. R.
2009
Groundwater monitoring networks
.
US geological survey
,
Nevada, USA
. In
Groundwater
(
Silveira
L.
Usunoff
E. J.
, eds). Volume
II
, pp.
441
.
Saaty
T. L.
1977
A scaling method for priorities in hierarchical structures
.
Journal of Mathematical Psychology
15
(
3
),
234
281
.
Saaty
T. L.
Alexander
J. M.
1989
Conflict Resolution: The Analytic Hierachy Approach
.
Praeger, RWS Publications
,
New York
.
Sharif
M.
Swamy
V. S. V.
2014
Development of LINGO-based optimisation model for multi-reservoir systems operation
.
International Journal of Hydrology Science and Technology
4
(
2
),
126
138
.
Song
J.
Yang
Y.
Chen
G.
Sun
X.
Lin
J.
Wu
J.
Wu
J.
2019
Surrogate assisted multi-objective robust optimization for groundwater monitoring network design
.
Journal of Hydrology
577
,
123994
.
Stefanakis
A. I.
2020
The fate of MTBE and BTEX in constructed wetlands
.
Applied Sciences
10
(
1
),
127
.
Vaezihir
A.
Zare
M.
Raeisi
E.
Molson
J.
Barker
J.
2012a
Field-scale modeling of benzene, toluene, ethylbenzene, and xylenes (BTEX) released from multiple source zones
.
Bioremediation Journal
16
(
3
),
156
176
.
Vaezihir
A.
Zare
M.
Vesali
Z.
2012b
Designing of Groundwater Monitoring Network Augmentation at Shiraz Oil Refinery
.
M.Sc. Thesis
,
Department of Earth Sciences, Shiraz University
,
Shiraz
,
Iran
.
Vanderbei
R. J.
2014
Linear Programming: Foundations and Extensions
, 4th edn.
Springer
,
Berlin
.
Yekom Consulting Engineers
2012
Semi-detailed studies of groundwater of plains under administration of regional water authority
.
East Azerbaijan Regional Water Authority
.
pp. 1–87 (in Persian)
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).