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

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.


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. ). 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. ). 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 ).
This goal can be achieved by establishing a groundwatermonitoring 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. ). 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 ). Groundwater-monitoring networks are constructed to meet regulatory requirements, to monitor groundwaterlevel fluctuation and groundwater quality variations through time, to determine contamination sources and to the evaluation of the efficiency of remediation (Hassan ).
Data, which can be collected from a monitoring program, are the direction of groundwater flow, thickness of the unsaturated zone, monitoring of short-and longterm groundwater-level fluctuations, tracking of pollution caused by released contaminants and their fate and transport (Rosen ).
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 selec- (Ohmer et al. ). 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 ), 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 ). can be used (Cunningham & Schrage ). 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 ). The design and optimization of groundwater-monitoring networks have widely been discussed in the literature over the last decades. Vaezihir et al. (a, b) 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. (), 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 () 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.
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 groundwatermonitoring 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.

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.
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 ). The general direction of groundwater flow in the area is from the southeast to the west and northwest (Figure 2).

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.  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 (w i ) and each node is covered (y i ) 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 ): x j ¼ (0, 1) for all j ∈ J P (4) y j ¼ (0, 1) for all i ∈ I P (5) where 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.
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. Algebraically, LP can be formulated in the following form (Vanderbei ): Subject to: X n j¼1 a ij X j b i i ∈ {1, : . . . :, m} X j ! 0 j ∈ {1, . . . , n} (Non À negativity constraints) (9) where x 1 …, x n 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.  (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 ).
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 () 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 (

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.

). 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 ).

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 Linear combination of all factors can be written as the following equation: where S is the total index, w i is the weight of factor i and x i 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.

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.

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: where y i , w i and x i 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.

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

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).
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). 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. ).
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.  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.

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.