ORIGINAL PAPER

Estimation of Curie point depths and heat flow from Ardebil

province, Iran, using aeromagnetic data

Allahyar Khojamli

& Faramarz Doulati Ardejani

& Ali Moradzadeh

Ali Nejati Kalate

& Amin Roshandel Kahoo

& Soheil Porkhial

Received: 12 March 2015 /Accepted: 26 February 2016 /Published online: 27 April 2016

The Author(s) 2016. This article is published with open access at Springerlink.com

Abstract This study attempt s to estimate the Curie point

depth (CPD) using the centroid and forward modelli ng of

the spectral peak methods in the Sabalan geothermal field in

Ardabil, NW of Iran. The reduced-to-pole aeromagnetic data

were divided into 18 overlapping blocks of the sizes of

100 × 100 km. In the centroid method, the average depth to

the top of the deepest crustal block, Z

, was first computed by

linear fitting to the second longest wavelength segment of the

power spectrum of aeromagnetic data. Then, depth to the cen-

troid of the deepest crustal block, Z

, was computed by linear

fitting to the longest wavelength segment of the power spec-

trum of the aeromagnetic data. The depth to the magnetic

bottom was obtained from Z

=2Z

− Z

. In the forward

modelling of the spectral peak method, the modelled spectra

fitted to the observed spectrum iteratively and Z

and Z

were

finally estimated. According to the obtained results, the CPD

varies from 10 to 18.6 km. The Curie temperature of magne-

tite was used to determine the thermal gradient and the heat

flow in the area. The study area is found to have a great energy

potential in the west, northwest a nd the southwest of the

Sabalan with shallow CPD, high geothermal gradient and heat

flow.

Keywords Curie point depth

Aeromagnetic

Geothermal

Sabalan

Spectral analysis

Introduction

The Curie point (approximately 580 °C for magnetite at atmo-

spheric pressure) is the temperature at which the spontaneous

magnetization vanishes and magnetic minera ls show paramag-

netic suscep tibility. The depth at which temperature reaches the

Curie point is assumed to be the bottom of the magnetized

bodies in the earth crust. Curie point temperature varies from

region to region depending on the geology and the mineralo g-

ical content of the roc ks. Therefore, one can normally ex pect

shallow Curie point depth (CPD) at the regions which have

geothermal potential, young volcanisms and thin crust (Aydin

and Oksum 2010). The assessment of the variations in the

Curie depth of an area can provide valuable information about

the regional temperature distribution at depth and the potential

of subsurface geothermal energy (T selentis 1991).

The idea of using aeromagnetic data to estimate CPD is not

new, and it has been widely applied to various parts of the

world. Bhattacharyya and Leu (1975) mapped Curie point

isothermal surface for geothermal reconnaissance of the

Yellowstone National Park in USA. In this area, CPD was

estimated 4–8 km. Tselentis (1991) calculated CPD in

Greece from aeromagnetic and heat flow data. Tselentis’sob-

jective was to understand the nature and extent of the regional

* Allahyar Khojamli

a.khojamli.g@gmail.com

Faramarz Doulati Ardejani

fdoulati@ut.ac.ir

Ali Moradzadeh

a_moradzadeh@ut.ac.ir

AliNejatiKalate

nejati@shahroodut.ac.ir

Faculty of Engineering, Malayer University, Hamadan, Iran

School of Mining, College of Engineering, University of Tehran,

Tehran, Iran

Faculty of Mining, Petroleum and Geophysics, Shahrood University

of Technology, Shahrood, Iran

Renewable Energy Organization of Iran (SUNA), Tehran, Iran

Arab J Geosci (2016) 9: 383

DOI 10.1007/s12517-016-2400-3

geothermal system at a depth beneath the area of Greece by

constructing the Curie isotherms. The results of his investiga-

tions revealed that the CPD varies considerably beneath

Greece, reaching 20 km towards western Greece and about

10 km beneath the Aegean. In East and Southeast Asia, CPD

was determined based on the spectral analysis of magnetic

anomaly data by Tanaka et al. (1999). In this study, they used

many heat flow data from the boreholes. The estimated CPD

for this area using centroid method varied from 9 to 46 km. In

addition, they predicted CPD from heat flow data. The CPD

estimated from the heat flow data were very similar to the

results of the CPD analysis of magnetic data. Dolmaz et al.

(2005) concluded that the study of earth crust’s thermal struc-

ture in SW of Turkey is useful to determine modes of defor-

mation, depths of brittle and ductile deformation zones and

regional heat flow variations. Karastathis et al. (2010)found

the deep origin of the geothermal fields and volcanic centres in

central Greece, by combining a trav el-time inversion of a

micro-seismic dataset together with a CPD analysis based on

the aeromagnetic data. They also found that a possible magma

chamber can be presumed by detecting a low seismic velocity

volume at depths below 8 km and the CPD estimation at about

7–8-km depth as well.

Bansal et al. (2011) estimated the bottom depth of magnetic

sources in Germany using aeromagnetic data. At first, they

proposed a modified centroid method to estimate the depth

to the bottom of magnetic sources. To assess the calculated

bottom depth of magnetic sources, the results were then com-

pared with the heat flow density data. Saleh et al. (2012)

estimated CPD and heat flow map for Northern Red Sea rift

of Egypt. Their aim was to map the CPD based on the spectral

analysis of the aeromagnetic data. The CPD varied from 5 to

20 km. The shallowest CPD of 5 km (associated with the high

heat flow) was suggested a promising area for geothermal

exploration. Eletta and Udensi (2012) investigated the CPD

isotherm from the aeromagnetic data to prepare a preliminary

potential map of geothermal resources in the Eastern Sector of

Central Nigeria. They showed that the high prospect areas are

located in the south-west parts of the study area. Obande et al.

(2014) applied spectral analysis of aeromagnetic data for geo-

thermal prospecting in the north-east Nigeria. They estimated

the top and the centroid depths of magnetic source from the

power spectrum. The obtained results were subsequently used

to estimate the bottom depth. The range of CPD varies from 6

to 12 km according to the heat flow and CPD values of the

study area wherein the highest heat flow value and the

shallowest CPD occurred near the thermal springs. The

Wikki warm spring area was found to have a great energy

potential with a shallow CPD and very high heat flow values.

The geological and geophysical evidences together with

the presence of several hot water springs in Ardebil province

in the NW of Iran indicate that the area could have a high

geothermal energy potential. Besides, the review of the

published materials shows that no comprehensive aeromag-

netic data analysis exists to prove the geothermal potential of

the region. So, any study regarding to locate geothermal po-

tential zones in such a vast region is highly important in the

early stage of a geothermal exploration program. Therefore,

this paper attempts to apply Centroid depth and forward

modelling of the spectral peak methods of the aeromagnetic

data to determine CPD in the main part of the Ardebil prov-

ince particularly around the Sabalan mountain area. The heat

flow values are then estimated and mapped to assess further

geothermal zones.

Geological settings

Ardabil province is a famous tourist destination in Iran. Its

pleasant climate especially in spring and summer seasons is

always worthy for most visitors and residents. Several hot

springs, with temperatures varied between 20 and 85 °C, exist

in Ardabil, which there are mostly around the Sabalan

Mountain (Mt. Sabalan). The Sabalan geothermal area

(Figs. 1 and 2) which is now under investigation for geother-

mal electric power generation lies at the NW of the Mt.

Sabalan (Ghaedrahmati et al. 2013). The area has been under

geoscientific exploration studies since 1978 (Fotouhi 1995).

Ardebil geology is diverse and complicated and has a long

evolution history. These features discriminate the area from

the other part of Iran. North of Ardabil is covered with older

alluvial, Clay, Marl and tuff intercalations. Surrounding region

around Mt. Sabalan is characterized by the predominance of

Quaternary terrace deposits (Dizu Formation); altered post-

caldera Pleistocene trachyandesitic domes, flows and lahars

(Kasra Formation); unaltered syn-caldera Pleistocene

trachydacite to trachyandesitic flows, domes and l ahars

(Toas Formation); and pre-caldera trachyandesitic lavas, tuffs

and pyroclastics (Valhazir Formation) (Fig. 2). The geologic

study of the Northwest of Sabalan confirmed that there are

two major types of structural setting: a set of linear faults

and several inferred faults (SKM, Sinclair Knight Merz

2005); the faults strike predominantly towards the northwest

and northeast (KML 1998). A northeast-southwest structural

trend is dominant in the south of Ardebil city. The main geo-

logical units exposed in this area include Miocene’saltered

tuff, tuff breccia, pumice, travertine, sandstone, shale, marl

and conglomerate and Eocene’s olivine basalt and

trachybasalt which overlay volcanic breccia and

trachyandesite of Eocene age.

However, the structural trend changes to northwest-

southeast direction in the further southern parts. This area

geologically contains a sedimentary sequence including

Cretaceous limestone, Jurassic’s shale and sandstone with in-

tercalation of dolomite which overlain by Eocene’svolcanic

383 Page 2 of 11 Arab J Geosci (2016) 9: 383

Fig. 1 Sabalan geothermal area

Arab J Geosci (2016) 9: 383 Page 3 of 1 1 383

breccia and trachyandesitic, trachybasaltic lava breccia and

lava flows of Quaternary age.

Data and methods

Aeromagnetic data of the area was obtained from the

Geological Survey of Iran. This data was corrected for the

International Geomagnetic Reference Field (IGRF 1976).

This data was collected in 1974–1975.

In this study, the centroid depth and forward modelling of

the spectral peak methods of the aeromagnetic data were used

to determine CPD.

Bhattacharyya and Leu (1977) presented a method for de-

termining the centroid of rectangular parallel piped sources,

which it had been used in their earlier study (Bhattacharyya

and Leu 1975) to investigate the Curie depths of Yellowstone

Park. This method was further developed by Okubo et al.

(1985) who combined and expanded the ideas of the methods

to the purpose of geothermal exploration. If the magnetization

of a set of two-dimensional (2D) bodies is completely random

and uncorrelated, the radial average of the power density spec-

tra of the total field anomaly, p(k), could be simplified as

follows (Blakely 1995; Stampolidis et al. 2005):

PkðÞ¼A

−2 kjjZ

1−e

− kjjZ

−ZtðÞ



ð1Þ

where A

is a constant and Z

and Z

represent the depths to the

top and bottom of the magnetic body, respectively. k denotes

the wave number of the magnetic field.

AccordingtoOkuboetal.(198 5), CPD (Z

)canbe

obtained in two steps. Firstly, the centroid depth (Z

)of

the deepest magnetic source is estimated from the slope

of the longest wavelength part of the spectrum divided by

the wave number using the following equation (Nwankwo

and Shehu 2015):

PkðÞ

¼ A− k

ð2Þ

where P(k) is the power density spectrum and A represents a

constant.

The depth to the top of the magnetic source is similarly

derived from the slope of high wave number portion of the

power spectrum as follows (Nwankwo and Shehu 2015):

ln PkðÞ



¼ B− k

ð3Þ

where B is a constant.

The depth to the magnetic bottom is then obtained from

=2Z

− Z

Fig. 2 Geologic map of the Sabalan geothermal field (modified from Ghaedrahmati et al. 2013)

383 Page 4 of 11 Arab J Geosci (2016) 9: 383

Forward modelling of the spectral peak method

Ravat (2004); Finn and Ravat (2004); Ross et al. (2004)and

Ravat et al. (2007) proposed forward modelling of the spectral

peak to better estimate the bottom depth by using Eq. 4.

PkðÞ¼Ce

− k

−e

− k



ð4Þ

where the constant C, non-depth-dependent term, can be ad-

justed to move the modelled curve up or down to fit the ob-

served peak. Location of the spectral peak and the slope in the

high wave number range are controlled by Z

and Z

,respec-

tively. The combination of both Z

and Z

controls the slope

immediately adjacent to the peak (Ravat et al. 2007). The

advantage of forward modelling is that it allows one to fit

iteratively the position and the width of the peak and match

the adjacent part of the slope more precisely and explore the

model space. Based on the fit of modelled spectra with the

observed, one may accept or reject the results more confident-

ly in this overall subjective process of fitting specific parts of

the spectra.

In both of these methods, CPD is computed in three steps:

(1) dividing the total magnetic field map into overlapping sub-

regions, (2) calculating the logarithm of power spectrum for

each region and (3) considering Z

=2Z

− Z

equation for

centroid method and calculating the basal depth or fit of the

modelled spectra with the observed spectrum and calculating

and Z

Having found the CPD, the heat flow values of the region

can be calculated using the following equation (Turcotte and

Schubert 1982):

q ¼ k

∂T

∂Z

ð5Þ

where q is the heat flow, k represents the coefficient of thermal

conductivity, and

∂T

∂Z

denotes the thermal gradient.

Data processing and analysis

Selecting the optimal dimensions of the sub-regions is very

important. The limited depth extent of the crustal magnetiza-

tion would be visible in magnetic maps, covering less than

100 × 100 km (Maus et al. 1997). Okubo et al. (1985, 2003)

suggested the optimal dimensions of the investigated square

window to be about ten times the actual target depth (Hisarli

et al. 2011). Connard et al. (1983) divided a magnetic data of

the Cascade Range, central Oregon, into overlapping blocks

(77 × 77 km) and calculated the radially average power spec-

trum for each block. Tanaka et al. (1999) divided the East and

Southeast Asia into sub-region data (approximately

200 × 200 km) and estimated the power density spectra for

each region. Blakely (1988) divided the area of Nevada into

blocks (120 × 120 km) in terms of magnetic or aeromagnetic

data and mapped the CDP of Nevada state.

Fig. 3 Selection of overlapping

blocks on the RTP map. Solid

circles indicate centres of blocks

Arab J Geosci (2016) 9: 383 Page 5 of 1 1 383

In the present study, the reduced-to-pole (RTP) aeromag-

netic data were divided into 18 overlapping blocks of sizes of

100 × 100 km (overlapped 50 % with the adjacent blocks)

(Fig. 3). It is widely acknowledged that the utilization of a

small window width may be a fundamental error in the appli-

cation of spectral methods for aeromagnetic interpre tation

(Nwankwo and Shehu 2015). If the source bodies have bases

deeper thanL/2π, they may not be appropriately resolved by

spectral method (Shuey et al. 1977). Therefore, data window

of 100 × 100 km possibly will satisfactorily resolve only depth

information to a depth of 15 km. This window size is based on

the fact that a computer code which has been written to deter-

mine a proper block size for calculating the radial power spec-

tra considers different window sizes varying from 50 to

400 km with a 10-km increasing step size. The absence of a

peak indicates that the peak lies at wave numbers lower than

the minimum resolved wave number; hence, a larger window

size is needed to compute the radial power spectrum and de-

tect the bottom of magnetic sources. The results show that a

peak is observed for some window sizes, but it shifts and

eventually disappears with increasing window size. The ap-

propriate block size of 100 × 100 km was then chosen so that

the maximum spectral peaks of the aeromagnetic data could

be visible in the power spectrum. However, few blocks with

CPD more than 13 km were recomputed with data window

increased to 200 × 200 km. The 2D power spectrum of

aeromagnetic data for each block was then computed using

the Oasis montaj software with fast Fourier transform

(FFT) method. Then, the effects of very deep regional

structures were removed using a first-order trend filter for

eachblock,andgridswereexpandedby10%usingthe

maximum entropy method to make the edges continuous.

The biggest advantage of 2D power spectrum is that the

depth of sources is easily determined by measuring the

slope of the power spectrum when the centroid method is

used (Saleh et al. 2012).

Results and discussion

This is the first time that spectral analysis of the aeromagnetic

data is being used to assess the thermal structure of the crust in

Ardebil province and mainly around Mt. Sabalan area. So, as a

representative example, the power spectrum plots for block

no. 17 in the north Sabalan are shown in Fig. 4. Once the

depth to the top of this block (Z

= 1.2 km) was estimated

using Eq. 3 and Fig. 4a, then Eq. 2 was applied to estimate

the centroid depth (Z

= 7.8 km) as shown by Fig. 4b. Finally,

CPD was calculated using Z

=2Z

− Z

= 14.4 km. Besides,

Fig. 4c shows an example of the spectral peak forward model-

ling method of the same block. In this figure, the calculated

power spectrum was fitted iteratively with the measured pow-

er spectrum using MATLAB software (Version 7.12.0.635,

R2011a). By using this method, the top depth of this block

= 4.1 km) and Curie depth (Z

= 14.1 km) were deter-

mined. The obtained results for the other blocks applying the

centroid depth and forward modelling of the spectral peak

methods are given by Tables 1 and 2,respectively.

Figure 5 also shows the CPD map of the study area for two

methods: (a) centroid method and (b) forward modelling of

the spectral peak method. In these figures, two shallow CPD

zones were observed. The first shallow CPD is identified in

the west and northwest of the Mt. Sabalan, and the second

zone is seen in the bottom left corner of the map which is

Fig. 4 Examples of power spectra for the block no. 17. a The depth of the

top Z

= 1.2 km is obtained by fitting a straight line through the second

longest wavelength spectral segment. b The depth to centroid Z

=7.8km

is obtained by fitting a straight line through the longest wavelength

portion of the spectra. The Curie point depth for this block is

=2Z

− Z

= 14.4 km. c The blue line represents the measured

power spectrum, and the red line is the result of forward modelling: its

depth to the top is 4.1 km, and bottom is 14.1 km

383 Page 6 of 11 Arab J Geosci (2016) 9: 383

located mainly in the East Azerbaijan province. In the centroid

method (Table 1), block no. 1 has the shallowest CPD (10 km)

and block no. 10 has the deepest CPD (16.9 km). Whereas, by

applying the forward modelling of spectral peak method, the

shallowest and deepest CPDs are assigned for blocks no. 1 and

no. 18, respectively (Table 2).

Table 2 Estimated CPD, geothermal gradients and heat flow by using the forward modelling of spectral peak method for the 18 blocks in the study

area

Block

number

Coordinates (UTM) Depth to top

(km)

CPD (km) Geothermal

gradient (°C/km)

Heat flow

(mW/m

)

Easting (m) Northing (m)

1 637,367.7 4,145,618 6.5 10.5 55.2 138.1

2 687,367.7 4,145,618 5.3 11.8 49.2 122.9

3 737,367.7 4,145,618 7.1 13.0 44.6 111.5

4 787,367.7 4,145,618 6.7 12.2 47.5 118.9

5 837,367.7 4,145,618 7.0 15.0 38.7 96.7

6 637,367.7 4,195,618 6.3 13.5 43.0 107.4

7 687,367.7 4,195,618 5.8 16.5 35.2 87.9

8 737,367.7 4,195,618 6.5 17.5 33.1 82.9

9 787,367.7 4,195,618 6.1 15.2 38.2 95.4

10 837,367.7 4,195,618 6.8 16.2 35.8 89.5

11 637,367.7 4,245,618 5.2 12.2 47.5 118.9

12 687,367.7 4,245,618 4.9 10.6 54.7 136.7

13 737,367.7 4,245,618 5.0 10.9 53.2 133.0

14 787,367.7 4,245,618 5.8 15.3 37.9 94.8

15 837,367.7 4,245,618 6.3 15.9 36.5 91.2

16 637,367.7 4,295,618 5.5 14.8 39.2 98.0

17 687,367.7 4,295,618 4.1 14.1 41.1 102.8

18 737,367.7 4,295,618 5.6 18.6 31.2 78.0

Table 1 Estimated CPD, geothermal gradients and heat flow by applying the centroid depth method for the 18 blocks in the study area

Block

number

Coordinates (UTM) Depth to centroid

(km)

Depthtotop

(km)

CPD (km) Geothermal

gradient (°C/km)

Heat flow

(mW/m

)

Easting (m) Northing (m)

1 637,367.7 4,145,618 6.3 2.6 10.0 58.0 145.0

2 687,367.7 4,145,618 7.4 2.5 12.3 47.2 117.9

3 737,367.7 4,145,618 7.2 2.4 12.0 48.3 120.8

4 787,367.7 4,145,618 7.8 2.7 12.9 45.0 112.4

5 837,367.7 4,145,618 9.1 1.7 16.5 35.2 87.9

6 637,367.7 4,195,618 7.9 2.4 13.4 43.3 108.2

7 687,367.7 4,195,618 8.8 2.6 15.0 38.7 96.7

8 737,367.7 4,195,618 9.1 1.9 16.3 35.6 89.0

9 787,367.7 4,195,618 8.9 2.2 15.6 37.2 92.9

10 837,367.7 4,195,618 9.5 2.1 16.9 34.3 85.5

11 637,367.7 4,245,618 7.8 2.2 13.4 43.3 108.2

12 687,367.7 4,245,618 6.9 2.4 11.4 50.9 127.2

13 737,367.7 4,245,618 6.8 2.3 11.3 51.3 128.3

14 787,367.7 4,245,618 8.2 1.8 14.6 39.7 99.3

15 837,367.7 4,245,618 8.1 1.6 14.6 39.7 99.3

16 637,367.7 4,295,618 8.9 1.9 15.9 36.5 91.2

17 687,367.7 4,295,618 7.8 1.2 14.4 40.3 100.7

18 737,367.7 4,295,618 9.2 2.2 16.2 35.8 89.5

Arab J Geosci (2016) 9: 383 Page 7 of 1 1 383

For determination of the surface heat flow, the value of

thermal conductivity (k)isgiven2.5Wm

−1

°C

−1

according

to the average thermal conductivity of the crustal rocks

(Stacey 1977). The geothermal grad ient of each block

was calculated by dividing 580 °C by the CPD. Heat flow

was then calculated by multiplying the geothermal gradient

by the thermal conductivity. Equation 5 was finally applied

to estimate the flow map of the study area (Fig. 6).

Figure 6aandTable1 illustrate the lowest heat flow value

(85.5 MW/m

) for block no. 10 and the highest heat flow

value (145 MW/m

) for block no. 1. In Fig. 6bandTable2,

one can see the highest heat flow value of 138.1 MW/m

occurring in block no. 1 and the lowest heat flow value of

78.0 MW/m

in block no. 18. Tables 1 and 2 give the

Fig. 5 The CPD maps of the

study area: a centroid method and

b forward modelling of the

spectral peak method (contour

interval is 0.5 km), thermal spring

(red solid circle), drilled well

(black square) and Mt. Sabalan

(yellow solid triangle)

383 Page 8 of 11 Arab J Geosci (2016) 9: 383

calculated CPD, heat flow and the geothermal gradient

values for a ll 18 blocks. The ma ps in Fig. 7 show the

geothermal gradient estimates each of these two methods.

It can be observed that a region with high temperature

gradient (Fig. 7a, 7b)andhighheatflow(Fig.6a, 6b) is

associated with the shallow CPD (Fig.5a, b). Results of the

both methods are rather the same. Also, the shallowest

CPD, h ighest hea t flow value an d maximum thermal

gradient are shown in the west, northwest and the south-

west of study ar ea. These results are confirmed by the

available well data, gravity and magnetotelluric measure-

ments. The CPD strongly varies according to the geologi-

cal conditions (Ross et al. 2006). The CPDs at volcanic and

geothermal areas are shallower than 10 km (Obande et al.

2014). Heat flow of about 80– 100 MW/m

indicates geo-

thermal anomalous conditions (Jessop et al. 1976). For the

Fig. 6 Heat flow maps in the

study area: a centroid method and

b forward modelling of the

spectral peak method (contour

interval is 5 mW/m

), thermal

spring (red solid circle), drilled

well (black square) and Mt.

Sabalan (yellow solid triangle)

Arab J Geosci (2016) 9: 383 Page 9 of 1 1 383

study area and area around Mt. Sabalan , the geothermal

anomalous conditions have been assigned to the most

blocks. As a result, t he study area is found to have a great

energy potential in the west, northwest and the southwest

of Mt. Sabalan with a shallow CPD, high geothermal gra-

dient and heat flow.

Conclusions

An attempt has been made to calculate the depth to bottom of

the magnetic so urces from the aeromagnetic data in the

Ardebil province in NW of Iran using spectral methods. The

CDP has been calculated by the centroid depth and forward

Fig. 7 Gradient maps in the

study area: a centroid method and

b forward modelling of the

spectral peak method (contour

interval is 2 °C/km), thermal

spring (red solid circle), drilled

well (black square) and Mt.

Sabalan (yellow solid triangle)

383 Page 10 of 11 Arab J Geosci (2016) 9: 383