Application of MCLP and LINGO methods to optimal design of groundwater monitoring network in an oil refinery site

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.

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.

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).

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 (NaN₃) 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.

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.

• w_i 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, J_p is the set of nodes j occupied by pre-existing facilities.

• N_i:⁠;

• d_ij is the shortest distance from node i to node j; S is the cover distance threshold (maximal service distance); x_j equals to 1 if a facility is located at site j; 0 otherwise; y_i 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).

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.

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).

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 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).

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).

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.

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.

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 groundwater-monitoring network was carried out in following steps.

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.

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

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.

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.

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.

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).

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).

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

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).

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).

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).

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).

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).

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.

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.

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.

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