Development of an updated geothermal reservoir conceptual model for NW Sabalan geothermal field, Iran

In this paper, a conceptual model has been developed for NW Sabalan geothermal field using exploration indicators. These indicators include the data and results of subsurface and surface investigations comprising geology, geophysics , hydro-geochemistry, hydrology and temperature and pressure distribution. The subsurface information was obtained from 10 deep exploration wells as well as from the results of previous studies in this field. All available data together with the stratigraphy and 1:20,000 geological map covering the study area were combined to produce a two-dimensional geological cross section. Also, a subsurface three-dimensional geological model was developed using available drilling logs. NW Sabalan geothermal field is one of the 18 detected potential areas in Iran that is located in Northwest of the country. This field includes a deep geothermal reservoir with a temperature range of 230–242 °C. The reservoir is covered by a cap rock with an approximate thickness of 500 m. There are four major geological units identified in the study area including Quaternary alluvium, fan and terrace deposits; Pleistocene post-caldera trachyandesitic Flows; Pleistocene syn-caldera trachydacitic to trachyandesitic domes and Pliocene pre-caldera trachyandesitic lavas, tuffs and pyroclastic. A hydraulic conductivity zone has been assessed by magnetotelluric surveys at deeper zones, suggesting that the main outflow direction is towards west and north of the area. The fluid chemistry is consistent with high chloride, neutral pH and immature liquids that partially equilibrated with host rock and they are classified as mixed water. The results of exploration and geological studies during a 10-year period have been integrated together to build a new overall conceptual model of the field.

fluid samples collected from wells, information on temperature and pressure conditions based on analysis of the available well logging data as well as other reservoir engineering information. A realistic conceptual model allows an explicit description of the main properties of a geothermal system such as hydraulic characteristics, groundwater flow, solute and heat transfer. A realistic and accurate conceptual model is probably not possible without a numerical model to prove its accuracy.
Northwest Sabalan (NW) geothermal field is one of the eighteen fields of identified potential and the fourth main high temperature areas in Iran that is situated in the northwest of Iran in Ardabil Province. Its distance from Tehran is 859 km (Fig. 1a) (Noorollahi et al. 2009;Porkhial et al. 2015;Kosari and Sattari 2015). This field has practically been exploring and developing recently. The area has been under geological investigation since 1978 (Fotouhi 1995).
The plan of electricity generation from NW Sabalan geothermal field has been presented since 1994. Afterwards emphasis has been put on this field as a top priority area. From 1998 to 2005, the detailed geo-based study was directed by the joint cooperation of renewable energy organization of Iran (SUNA) and Sinclair Knight Merz Ltd (SKM) of New Zealand. The NW Sabalan geothermal field was finally distinguished as an important potential for power generation purpose (Noorollahi et al. 2009).
Accordingly, during the time period of 2002-2004, three deep exploration and delineation wells (NWS-1, NWS-3 and NWS-4) have been drilled. Besides, a shallow injection well named NWS-2 (Fig. 1b, PADs A, B and C) was drilled, in order to evaluate the subsurface geology and provide data for assessment and modelling the reservoir (Najafi and Ghobadian 2011; Porkhial et al. 2015). The location of these wells was determined based on the results obtained from a Magnetotelluric (MT) survey conducted in 1998 (SKM 2005a).
The wells NWS-1 and NWS-4 were the first deep wells drilled in Sabalan volcanic area (Porkhial et al. 2015;Kosari and Sattari 2015). A maximum temperature of 242 °C was measured in well NWS-1 at a minimum depth of 832 m. Subsequently, from 2008 to 2011, SUNA performed new geophysical exploration studies. Based on the new results, the wells NWS-5 (T max = 240 °C, D min = 1635 m) and NWS-11 (T max = 174.5 °C, D min = 2701 m) on Pad A and C, respectively, NWS-6 (T max = 238 °C, D min = 1144 m), NWS-7 (T max = 237 °C, D min = 1341 m) and NWS-10 (maximum temperature 242 °C, D min = 1890 m) on Pad D and NWS-8 and NWS-9 (maximum temperature 230.5 °C at a minimum depth of 2301 m), respectively, on Pad E were drilled (Fig. 1b).
A few studies have been carried out in the past to present a conceptual model for the NW Sabalan geothermal field. The conceptual/geological models presented by KML (1999b), Noorollahi and Itoi (2011a, b), Porkhial et al. (2010a, b), Khojamli (2012) and Sabzemeidani (2016) are all noteworthy. KML (1999a) built a hydrogeological geothermal model for the Mt. Sabalan volcanic complex based on multi-disciplinary combination of geological, geochemical and geophysical data. The heat source, exploration wells, meteoric waters, deeply circulated meteoric waters, magmatic volatiles and condensates are illustrated. Based on this model, one may explain that a large residual magma mass is shown to underlie the Sabalan caldera complex with intrusive apophyses developing upwards from the magma mass to shallow depths controlled by caldera ring-faults. KML (1999a) further explained that the thermal features might all have a common origin within the Sabalan caldera. Noorollahi and Itoi (2011a, b) have presented a conceptual model from the available three deep exploration wells data as a basis for their reservoir numerical model. The simulation predicted a power plant equivalent to be designed with a capacity of 50 MWe and an age of more than 30 years. Now construction of the first pilot plant is running with a capacity of 5 MWe, the produced fluid temperature of 86 °C and maximum flow rate of 58 l/s (average flow rate: 46 l/s) in the NW Sabalan Site. This plant is located between Pad A and Pad D with geographical coordinates x = 739,108, y = 4,238,580 and elevation of 2635 masl (Fig. 1b).
Literature review indicates that despite several researches related to the development of conceptual models for NW Sabalan geothermal field, a comprehensive conceptual model based on the combination of geological, geophysical, hydrogeochemical, hydrological and deep well data has not been yet presented to provide a basis for developing a numerical model for the study geothermal reservoir. The weakness of these works (KML 1999a, b;Porkhial et al. 2010a, b;Noorollahi and Itoi 2011a, b;Khojamli 2012;Sabzemeidani 2016) is that they only employed the information obtained from three exploration wells for the development of the conceptual model. Besides, they ignored some of important exploration indicators to be considered. Therefore, these models cannot fully explain the processes of Sabalan geothermal system, such as conduction and convection zones. In addition, they presented an incomplete explanation about the fluid flow cycle. Moreover, they did not define an exact boundary for the geothermal reservoir due to the lack of adequate data. In the present work, attention has been focused on the development of a conceptual model which takes all geothermal system variables into account and provides a visual representation of the geothermal system. The structure of the conceptual model presented here was completed by the use of geological and thermal data obtained from the new deep exploration wells. Conduction and convection zones and the top of the reservoir were added to this model by the incorporation of the geophysical results together with subsurface temperature distribution. To derive the conceptual model, Datamine Studio3 Version 3.21.7164 and ArcGIS Desktop software were used as tools.

Geology
NW Sabalan geothermal field is located in the northwest part of Mt. Sabalan. A geological map at a scale of 1:20,000 of Sabalan region is shown in Fig. 1b. Mt. Sabalan is a large stratovolcano comprising a wide-ranging central structure which is constructed on a likely tectonic horst of underlying intrusive and effusive volcanic rocks. According to KML (1999a), the output magma had a great role in forming a collapsed caldera about 12 km in diameter and a depression of about 400 m.
The Mt. Sabalan region sits on the south Caspian plate, which underthrusts the Eurasian Plate to the north and it is in turn underthrust by the Iranian plate, which causes compression in a NW direction. Besides, a dextral rotational motion resulted by the movement of Arabian plate beneath the Iranian plate makes the geological structures in the region more complicated. As a result of this tectonic structure, the Cenozoic geological history and stratigraphy of the region under consideration are complicated and the geological units in this region might have different structural characteristics (Mckenzie 1972;Amidi 1978;TBCE 1979;ENEL 1983;Nejad 1987;Manouchehri 1989;Emami 1994;Noorollahi et al. 2008).
Two major geological structures have been identified in the study area that comprise a set of linear faults and several ring-faults ( Fig. 1b) (KML 1999a;SKM 2005a, b;Noorollahi et al. 2008;Ghaedrahmati et al. 2013). The results of interpretation of Centre National d'Etudes Spatiales (SPOT) satellite images and aerial photographs showed a WNW trending structural zone passing through the Moeil Valley. The faults strike mostly towards the NW and NE (KML 1999a).

Stratigraphy
By increasing data due to adding new exploration wells and constructing Pads D and E, the realisation of the deep volcanic succession in the NW Sabalan geothermal field became easier. These data have been used in this paper to understand the stratigraphy in details. It is noted that the deepest geothermal wells in NW Sabalan are NWS-1 and NWS-3 with approximate depths of 3177 and 3197, respectively.
A three-dimensional geological model ( Fig. 2) was generated using both the available wellbore and surface geological data. Datamine Studio3 Version 3.21.764 software was used as a tool. Figure 3 shows a two-dimensional cross section including surface and subsurface geological data. The geological map (Fig. 1b), wellbore data and the threedimensional model (Fig. 2) were employed to build this section.
Monzonite basement rock was revealed at the depth of 1021 m of NWS-1 well in southeast. This rock was further approved at the depth of 2140 m of NWS-6. However, it appeared as the monzonite dykes at the depth zone of 1978-2700 m of NWS-9 well. The outcrop of monzonite extends as a bound, 5-10 km wide, for 100 km with a general NW strike. Mineralogical study of this rock showed minerals comprising plagioclase, orthoclase and slight amounts of quartz, biotite and minor hornblendes. Also, the wells NWS-4 and NWS-5 reached to Diorite at an approximate depth of 1840 m. This rock type appeared as old and young diorite intrusive in NWS-1.

Epa Formation
The thickness of Epa formation varies from a few hundred metres to 2200 m in the study area and it includes trachyandesite, trachybasalt, andesite, altered andesite and lahar. Name of this unit has been taken from the local geological map at a scale of 1:100,000 (Emami 1994). This formation has a separate weathering horizon at its upper surface. This unit identifies an unconformity with the overlying Valhazir formation (Fig. 3)

Valhazir formation
This formation has been introduced by Bogie et al. in (2000). Its thickness varies between 300 and 2000 m with a variable alteration zone (from 420 to 970 m). This formation includes Pliocence pre-caldera trachy-andesitic lava flows, tuffs and pyroclastic breccia as the oldest volcanic rocks identified in the study area. According to Bogie et al. (2000) these volcanic units are more fractured and affected by faulting than the younger units.
The subsurface studies identify that pyroclastics is dominant in the upper parts of the unit with lavas in lower portions (SKM 2005a, b).
The Valhazir formation is equivalent with Qpad of the geological map of Meshginshahr (Emami 1994). Samples of Valhazir formation taken through the exploration well is totally altered to clay, silica and pyrite with slight amounts of chlorite and Fe-oxides. This means that the initial lithology of Valhazir is not clear, hence a term of altered volcanic is used for this formation (Eshaghpour 2005).
Coarser-grained pyroclastics comprise lapilli tuff and tuff breccia. However, their fine cutting size cannot be easily recognised. Lithic clasts are dominant and include a variety of andesites with either biotite or hornblende as the main mafic minerals. The andesite at the base of the formation includes euhedral phenocrysts of plagioclase and hornblende in a dark and very fine-grained matrix (SKM 2005a, b).

Plat and Ngab formation
These units have been named by Emami (1994) and were identified in NWS-3 well, only. The main indication for the presence of these units in NWS-3 is the existence of olivine phenocrysts in the altered rocks with trachybasalts. Trachybasalts have been first reported in the stratigraphic sequence in the Ngab unit, which is indicated in the cross section of geological map presented by Emami (1994). The Ngab unit is underlined by the Plat unit (Fig. 3). Plat formation includes Pliocene Lahar, Crystalitic Tuffs and Andesitic Breccia with a thickness of 114 m (SKM 2005c).

Dizu formation
The Dizu formation has been identified to be the upper stratigraphic sequence at NW Sabalan geothermal field (Bogie et al. 2000). It is equivalent to the Q t1 and Q t2 units (Amini 1988).
This formation comprises terrace deposits displaying a mixture of poorly sorted sand, granules, pebbles, cobbles and boulders with some sand intercalations. The terrace deposits show a blend of fluviatile processes, mainly during flash floods and by mass flow from the surrounding slopes (SKM 2005a; Mousavi and Darvishzadeh 2010).
The rounded clasts of andesite, trachyandesite and subordinate trachydacite are dominant in this formation. The andesite includes phenocrysts of plagioclase and hornblende in a grey matrix. The trachyandesite also comprises plagioclase and hornblende phenocrysts, but further includes K-feldspar phenocrysts showing a brown matrix. The trachydacite consists of coarse phenocrysts of sanidine along with plagioclase and hornblende and displays a light grey vesicular matrix. The coarser grains are very poorly sorted and clasts range from granule to boulder size, with a sand matrix. They are uncemented and show no indication of hydrothermal alteration. However, the rare clasts of chalcedonic quartz are the erosional products of the hydrothermally altered sequences of the Valhazir formation (SKM 2005a, b).

Geophysics
Geophysical surveys including Direct Current (DC) resistivity, Transient Electromagnetic (TEM) and MagnetoTelluric (MT) were conducted in the Mt. Sabalan region by the Renewable Energy Organization of Iran (SUNA) in 1998. The surveys were carried out at 212 resistivity stations over an area of about 860 km 2 . 110 of these stations were located in NW of Mt. Sabalan. These surveys resulted in five anomalies that the largest one is located at Moeil Valley (Fig. 1b) (KML 1999a, b;Bromley et al. 2000;SKM 2005a;Noorollahi et al. 2008;Ghaedrahmati et al. 2013).
The resistivity of geological layers at shallow depths was provided by DC resistivity data, and the magnetotelluric soundings with a spacing of 50 m were employed to provide resistivity data for deeper layers. In order to provide resistivity data for those layers at a depth zone of 100-300 m and for the treatment of "static-shifts effects", the TEM surveys have been conducted. The natural-source MT method with a frequency range of 8 kHz-0.02 Hz was employed to estimate the exploratory drilling depth. The MT method identified an area in the Moeil valley (NW Mt. Sabalan) with a low resistivity anomaly which is accessible at a depth range of 300-3000 m (Bromley et al. 2000;Noorollahi et al. 2008). Three deep exploratory wells (namely NWS-1, NWS-3 and NWS-4) were consequently drilled between 2002 and 2004 (at Pads A, B and C, Fig. 1b). Their locations were chosen based on a 4 Ωm low resistivity anomaly zone. However, this drilling program was not successful to identify the exact location of the upflow zone (SKM 2005a).
Subsequent interpretation of good quality MT data by SKM (2003) and Talebi et al. (2005) revealed a new location for outflow zone near the well NWS-3. The anomalous zone with low resistivity extended to southeast of NWS-1 (Pad D, E, Fig. 1b) that suggests a location for likely upflow zone. They recommended further investigation of the upflow zone location applying low-frequency magnetotelluric data and further drilling operation in this location.
A new MT study incorporating 78 deep soundings with a frequency range of 385-0.0005 Hz was conducted by EDC (EDC 2008) in 2007, to determine the location of likely reservoir and possible drilling targets for the development of a geothermal power plant. The results showed a shallow resistivity anomaly as the likely centre of geothermal system on the east of Pad E and west of the young trachyandesite domes of post-caldera volcanic formation. Due to inadequate coverage of MT stations, especially, at the east of the study area, another MT survey was designed and conducted in 2009 for more identification of the geothermal system (EDC 2010). Based on the joint interpretation of new MT sections and two-dimensional geological models (EDC 2010), a final analysis was performed between years 2007 and 2009 to confirm the presence of a hotter zone in the east side of Pad D. Ghaedrahmati et al. (2013) developed a three-dimensional invention code using a long-period MT data to obtain a realistic resistivity model. The resistivity slices of the different horizons from their model are shown in Fig. 4. A low resistivity zone [zone (1), on the resistivity slice map related to the depths below 300 m] has extended from northwest to the centre of the area. This zone is almost parallel to the general structural trend and faults of the study area. This zone has been identified as an alteration zone due to the upward thermal flows along the faults and permeable structures. These structures have been delineated from wells NWS-1 and NWS-4 and further validated by the presence of thermal springs (Fig. 1b).
The resistivity cross section given by EDC (2010) with a direction almost similar to that of geological section (Fig. 1b) was considerably modified and is shown in Fig. 5. The indicators of alteration zone identified in wells of Pad D and Pad E are also illustrated in this figure. Joint interpretation of Figs. 4 and 5 reveals that the geothermal resource is associated with the hydraulic conductivity zones (3) and (4) at an approximate depth of 4600 m. This deep conductive zone could be related to fluid flow within the zones consisting of high fracture density.
Also, two additional hydraulic conductivity zones can be highlighted in southern and south-eastern area in the 1000 m depth resistivity slice map [zones (5) and (6) in Fig. 4, first row, right figure]. The results of 3-D modelling approximated a depth more than 1500 m for these two conductive zones (see Fig. 4, second row, left figure). These zones Fig. 4 The resistivity slices at different depths (modified from Ghaedrahmati et al. 2013), the black dots on maps show MT stations and the labels 1, 2, 3, 4, 5 and 6 display conductive (anomaly) zones related to geothermal resource obtained from 3-D modelling results. New location of the hotter zone (as an ellipse) and its previously identified location by EDC (2010) (as a circle) are also illustrated on the 4600 m depth slice map are probably altered due to interaction of thermal fluids and country rocks along the main fractured zones and faults (EDC 2008(EDC , 2010Ghaedrahmati et al. 2013).
Hydraulic conductivity zone is one of the main indicators of hydrothermal activity in Mt. Sabalan area with resistivity values less than 20 Ωm which have been identified at elevation zones varying from 1200-2800 masl. The presence of hydraulic conductivity zones at deeper horizons suggests that the main outflow directions in Mt. Sabalan area are towards the west and north. These conductive zones have about 600-1000 m thickness below Pad D (Fig. 5).
In addition, the data obtained from well NWS-7 show that smectite and illite-smectite units are coincident with the hydraulic conductivity zones and also epidote unit is coincident with the increase in resistivity values (>30 Ωm). From the well NWS-8 data, it is obvious that smectite, chlorite and illite units have resistivity values less than 30 Ωm, whereas illite and smectite units are coincident with the resistivity contour lines above 30 Ωm (modified from Ghaedrahmati et al. 2013).

Geochemistry
The surface geochemistry of Mt. Sabalan area was first time interpreted by Stefansson (1989) and further reported by Fotouhi (1995). Subsequently, the geochemistry of hot springs has been studied by Khosrawi (1996). Similar investigations have been carried out by Stelbitskaya and Radmehr (2010) and Kosari (2011).
The natural thermal features at NW Sabalan have been discovered as hot springs, most of which are located on the outskirts of the Moeil valley with a temperature range of between 25 and 85 °C. This area is characterised with the most thermal features in Mt. Sabalan (Fig. 1b).
The spring waters are acid SO 4 , acid Cl-SO 4 and natural Cl-SO 4 types (Table 1). Hydrothermal alteration can be seen at the land surface. The main hydrothermal alteration is associated with the Valhazir formation (Rahmani 2007).

Springs chemistry
All thermal features in NW Sabalan arise mainly from the gravels of the Dizu formation (Fig. 1b). There is no spring reported in downstream areas at lower elevations. Their water is characterised with a simple dilution trend indicating mixing with varying amounts of shallow groundwater (Table 1). They show a strong seasonal cyclic variation in flow rate; however, very slight seasonal deviation in temperature or chemistry can be seen. Despite the elevated Cl concentration, isotopic data categorised the waters to lie on the local meteoric water line. Mg concentrations exhibit a weak inverse correlation with Cl. This may suggest that geothermal fluids containing elevated Cl concentrations are diluting with the shallow groundwater aquifer consisting of high Mg concentrations.
Assuming the fact that all of the springs are not boiling, quartz geothermometer without stream loss (Fournier 1977;Fournier and Potter 1982) may be acceptable and a temperature range of 80-156 °C was calculated using this geothermometer. The Na-K geothermometer resulted in temperatures greater than 250 °C. However, elevated Ca and Mg concentrations reject application of this geothermometer. The Na-K-Ca geothermometer calculated a temperature range of 75.5-213 °C. Accordingly, recent temperature estimation applying various geothermometers by the authors is given in Table 2. As well shown, the reservoir temperature at NW Sabalan varies between 88 and 244 °C.
According to the chemistry of the springs water and the different temperature ranges estimated by various geothermometers and also comparison them with the measured temperature of the reservoir (242 °C, in well NWS-10), it seems that the agreement between the calculated and measured reservoir temperature is somewhat close.

Reservoir chemistry
Although not given, several samples were collected from well NWS-1 in May and June 2004 and from well NWS-4 in September 2004. The samples were collected by SUNA. The samples were both untreated and acidified water. Chemical analyses of pH, Cl, HCO 3 , SO 4 , Ca, B and CO 2 were carried out in the site laboratory of SUNA and for Li, Na, K, Ca, Mg and SiO 2 in laboratory of GNS Wairakei (New Zealand). Stable isotope analysis (δ 18 O and δ 2 H) was conducted in the laboratory of GNS Wellington (Stelbitskaya and Radmehr 2010).
In the second stage of development and reservoir assessment, the discharged fluids from four wells, including NWS-6D, NWS-7D, NWS-9D and NWS-10D were analysed. The water samples have been classified as high chloride, neutral pH and mature liquids showing partial equilibrium with host rock. The analytical data of liquid samples comprise Ca, Na, K, Mg, Li, Fe, Mg, Mn, B, Cl, F, SiO 2 , SO 4 , CO 2 , H 2 S, HCO 3 and NH 3 concentrations in the liquid phase of Webre and Weir box samples. The steam samples were analysed for all gas components containing CO 2 , H 2 S, He, H 2 , Ar, O 2 , N 2 , CH 4 and NH 3 (Kosari and Sattari 2015).

Liquid chemistry
Based on SKM reports on production wells NWS-1 (SKM 2005a) and NWS-4 (SKM 2005b) and further interpretation by Rahmani (2007), the discharging fluid chemistry of wells NWS-1 and NWS-4 was classified as an alkaline-Ph, medium salinity and sodium chloride water. The chemistry of discharging fluid was almost identical for most of the chemical concentrations except Ca, which slightly varies from well NWS-4 (25-30 ppm) to well NWS-1 (15-30 ppm). In well NWS-1, Ca concentration decreased compared with Cl over the first week of discharging. Kosari (2011) suggested that the initial water discharged from the wells may have originated from deeper and different inflow zones with relatively different compositions apparently with higher Ca concentration.
The calculated calcite saturation of deep liquid phase in the aquifer water through wells NWS-1, NWS-4, NWS-6, NWS-7, NWS-9 and NWS-10 (Fig. 6a, b) revealed that the deep liquid is generally super-saturated at lower temperatures than reference temperatures and under-saturated at temperatures more than reference temperatures (see reference temperatures in Table 3) (SKM 2004a;Kosari and Sattari 2015). Figure 7 shows Na-K-Mg ternary diagram (Giggenbach 1991) related to the springs of NW Sabalan. It is well indicated that most of the samples are plotted in the partial equilibration and mixed waters domains. Applying two geothermometers including K-Mg and Na-K (Giggenbach 1988) on this geo-indicator reveals that the deep fluid temperatures range is from 235 to 320 °C and between 285 and 335 °C, respectively. Table 3 shows the reservoir temperatures calculated by several geothermometers and compared with the reference temperatures (Strelbitskaya and Radmehr 2010;Kosari and Sattari 2015). This table further indicates that there is a close agreement between temperatures calculated by quartz geothermometers and the reference temperatures.

Gas
The chemistry of discharging stream from the wells indicated that about 98% of the total gas content is CO 2 , while H 2 S, N 2 , Ar, H 2 and CH 4 form the remaining steam components (Kosari 2011). The temperature calculated by gas geothermometers, particularly H 2 S geothermometer changes from 269.1 to 277.1 °C (with an average of 274.2 °C), is fairly well correlated with that of Na-K geothermometer (279 °C). However, CO 2 geothermometer estimated temperatures ranging from 287.3 to 308.9 °C (having an average of 301.8 °C) which is higher than those calculated applying all the other gas and solute thermometers ( Table 3). The elevated CO 2 concentration in deep liquid is the main reason for this overestimation (Kosari and Sattari 2015).

Hydrology
In NW Sabalan area, the magmatic intrusions are permeable due to the abundance of fractures, typically in the upper parts. Field investigations have shown that the volcanics of the Oligo-Quaternary succession made up mainly of lava at shallow depths are also permeable. Pyroclastic deposits and altered zones have medium-to-low permeability (Khosrawi 1996). The permeability of these rock units eases the movement of recharging water to the subsurface. However, the hydrogeology of the NW Sabalan is mostly controlled by faults.
The average annual atmospheric temperature is about 8 °C while in higher elevations, it decreases to 0 °C. The average temperature is very low from December to February and it reaches to 0 °C or even below (−15 °C in January). The maximum value of average monthly temperatures is obtained in May and April. The average monthly precipitation has been recorded in synoptic stations between 1996 and 2015. The precipitation occurs mainly in the form of snow in the winter. Based on available data, the maximum average monthly precipitation is 72.79 mm in May and its minimum is 12.14 mm in August (Fig. 8).

Subsurface distribution of temperature and pressure
The reservoir thermodynamic data including temperature and pressure were measured from the exploration boreholes. In a geothermal system, there is a need to specify the temperature distribution in assessing reservoir characteristic. To achieve the goal, it is necessary to measure downhole temperature under an equilibrium condition and to be described under the stable state of reservoir. To study the reservoir temperature changes vertically and horizontally, the contour plots and vertical cross sections have to be provided, accurately. These plots are useful to show how the interaction of hot and cold fluids occurs, so they are very important in developing a proper framework for the hydrogeological model of the system (Abdollahzadeh Bina 2009;Porkhial et al. 2015). Figure 9 shows the stable temperature profiles for ten wells including NWS-1, NWS-3, NWS-4, NWS-5, NWS-6, NWS-7, NWS-8, NWS-9, NWS-10 and NWS-11. According to a measured water level equal to 2413 masl at well NWS-6, a representative 'Boiling-Point-for-Depth' (BPD) curve is also shown. Since, all of the temperature profiles are below the BPD curve, one may say that the reservoir in the study area contains water only. At an altitude of about 1800 masl, behind the casing of well NWS-1, the condition is near the BPD curve. It means that there is probably a two-phase fluid condition. For wells NWS-6, NWS-7 and NWS-10, the temperatures below the altitude of +1500 masl are near the BPD curve. The temperature reverse mode is slightly observed between +700 and −200 masl. Below the altitude of −200 masl at well NWS-1, +600 masl at well NWS-3 and +300 masl at well NWS-11, the temperature increases with depth. The highest temperature can be estimated at about 240-243 °C in well NWS-1 and in the altitude of about +1800 masl and also at well NWS-10 at an altitude of about 850 masl. There is a deep upflow zone between the wells NWS-1, NWS-6, NWS-7 and NWS-10, vertically and in southern direction.
By the use of measured temperatures in all exploration wells, temperature contours were plotted in a cross section along the section AB of Fig. 1b. The result is shown in Fig. 10a. As illustrated in this figure, the temperature increases from north to south. This means that there is possibly a high temperature upflow zone in the southern part of the area under investigation which is located below the Pad D. At a depth zone of about 1300-1750 masl, between wells NWS-6, NWS-7 and NWS-10 and NWS-1 located between Pad D and Pad A, the highest temperature of the reservoir (>240 °C) can be estimated. Also, a hot outflow from there originates and flows to SE-NW direction towards the well NWS-3. Besides, the cap rock with a thickness of about 500 m can be well recognised in the range of about 2300-1800 masl. The main flow below the drilled area covers at least a 2200 m thick convection zone approximately below the altitude 1700 masl.
Several faults, namely NW3, NW5, NE5, NNW2 and NW4 build the upflow zone for the reservoir system. Figure 10b shows the cross section of estimated pressure using available data obtained from 10 deep wells. This cross section was built using the pressure profiles of the wells which are equivalent to the pressure formed at the major feed zones, only. Figure 11 shows the final conceptual model for the cross section AB of Fig. 1b. This model has been developed based on the geological, geophysical, geochemical, temperature data described in the previous sections. The results of geophysical survey revealed that the geothermal reserve is associated with the hydraulic conductivity zones below the Pads E and D, 1000 m below sea level. In addition, the hydraulic conductivity zones were identified below the Pad D and E by resistivity values of less than 20 Ωm at elevation zones varying from 1200 to 2800 masl. The results of the temperature profiles near the reservoir show that the top of the reservoir is at about 1800 masl, below which there are little temperature changes in the convecting reservoir. The model includes a deep geothermal reservoir with a temperature range of 230-243 °C. Faults NW3 and NE5 possibly control the conductivity of the reservoir. Moreover, they have important role in the creation of upflow zone below the wells NWS-6, 7, 10 (Pad D) where the heat and pressure differences cause a single-phase flow travel through these faults and the other fractures by convection phenomenon. It is noteworthy to say that there is a good Fig. 10 Distribution of subsurface pressure and temperature for cross section AB of Fig. 1b (a temperature, b  pressure) correlation between the presence of hot springs and faults in the study area. It seems that the upward flow takes place through faults NNW2 and NE5 and NE1 and eventually discharges to surface as hot springs. According to the results of stable temperature profiles (Fig. 9) and distribution of subsurface temperature (Fig. 10a), it is remarkable to express the idea that the infiltration of meteoric water creates a shallow cold zone and this mixes with the upflow zone moving upwards to the shallower zones in subsurface. By penetrating meteoric waters to deeper zones, the geothermal system is fed and the fluid flow circle is eventually completed.

Conclusions
Identifying physical and chemical processes involved in a particular problem and creating an accurate conceptual model are the most critical issues in the development of a successful numerical model. An appropriate conceptual model can be used to describe the crucial features of geological situations and define the principal processes in a geothermal system. In this study, a conceptual model has been presented based on the available geological, exploration and drilling data for NW Sabalan geothermal field. The model includes a deep geothermal reservoir with a temperature range of 230-243 °C. This model can provide useful information regarding temperature, pressure, outflow zone and fluid flow circulation within the system that can be successfully used to improve the previous numerical models and provide a basis for sustainable development of geothermal field for imminent targets. However, timely updates of this conceptual model based on new surface and subsurface exploration data are critical for successful development planning, well siting and resource assessment of NW Sabalan geothermal field. Faults NW3 and NE5 possibly control the conductivity of the reservoir rocks. Furthermore, they have important role in the creation of upflow zone below the wells NWS-6, 7, 10 on pad D.