Electronic Journal of Polish Agricultural Universities (EJPAU) founded by all Polish Agriculture Universities presents original papers and review articles relevant to all aspects of agricultural sciences. It is target for persons working both in science and industry,regulatory agencies or teaching in agricultural sector. Covered by IFIS Publishing (Food Science and Technology Abstracts), ELSEVIER Science - Food Science and Technology Program, CAS USA (Chemical Abstracts), CABI Publishing UK and ALPSP (Association of Learned and Professional Society Publisher - full membership). Presented in the Master List of Thomson ISI.
Volume 22
Issue 2
Geodesy and Cartography
DOI:10.30825/5.ejpau.171.2019.22.2 , EJPAU 22(2), #01.
Available Online: http://www.ejpau.media.pl/volume22/issue2/art-01.html


H.S. Kutoglu1, M. Berber2, H. Kemaldere1
1 Geomatics Engineering Department, Faculty of Engineering, Bulent Ecevit University, Zonguldak, Turkey
2 Department of Civil and Geomatics Engineering, California State University, Fresno, CA, USA



The source of deformations in one of the most densely populated residential areas of the city of Zonguldak is studied using InSAR techniques. During the investigations, it is discovered there had been mining acitivites underway under the city during the time period that the InSAR images were taken. Approximately 8 years of this time period is examined, and in total a maximum 41 cm subsidence is discovered in the area of coal mining. In addition, the vertical motion coefficient for the area of interest is calculated as 1.6.

Key words: Geodesy, remote sensing, InSAR, subsidence.


Zonguldak is located at the Western Black Sea coast of Turkey, 360 km east of Istanbul and 270 km north of Ankara (Fig. 1). The topography around the city is very undulating and steep; 56% of the city’s land is surrounded by mountains, and only 29% where the slope is less than 20% and suitable for settlement and agriculture. Also, very heavy forests cover the land surface in the immediate vicinity of the city centre (Zonguldak Province Environment and Forestry Directorate, 2006).

Fig. 1. Location of Zonguldak city and hard coal basin (source: Google Earth)

The city, with a population of 200,000, relies on many sources of income, including forestry and mining. However, it is the main industry of coal mining that makes Zonguldak the most famous mining region of Turkey. Underground coal extraction in the basin was initiated in 1848. Currently, it is the Turkish Hard Coal Enterprise (TTK) company authorized for mine production in the basin. According to the official records of TTK, the hard coal production is about 2.5 million tons per year, and since 1848 has produced a total of 400 million tons. There are widespread coal seams located between the levels of -155 m and -550 m under the city (Turkish Hard Coal Enterprise, 2008).

Due to the above-mentioned mining activities, major subsidence events have occurred at some locations of the city in the past [3, 4]. Some locations are still suffering from subsidence phenomena, and some locations are at risk of possible damage because of subsidence. For the safety of human life and real property, monitoring the temporal evolution of the subsidence effects for the city of Zonguldak is very crucial. Also, such a monitoring program can facilitate defining the locations at risk, offer early warning before major occurrences, and guide proper urban planning for the future. For this purpose, a research project was initiated in 2005 to monitor the earth surface in and around the city against subsidence events caused by mining activities. Until the end of 2006, analyses were performed on RADARSAT data whereby a great portion of vertical movement was discovered in some parts of Bahcelievler district, which is a commonly favored and heavily inhabited location of the city. Whereas the data of the 2006 analyses found extensive movement, in the analyses conducted using JERS data in 1995 there was not any apparent movement in this location. When the source of the motion was investigated, it is found out that in 2005 mining activities began just underneath this place. These coal mining activities continued until 2011 (Fig. 2).

Fig. 2. Mining activity under the area of interest. The solid areas in red shows the driven coal panels, and the colored lines represent the mine seams at different levels changing between -550 and -230 m (source: Google Earth)

In order to thoroughly monitor vertical movement stages caused by coal production, analyses have continued using PALSAR, ENVISAT and TerraSAR-X images [1, 4, 5]. In this study, focus is on Bahcelievler district under which coal production had continued between 2005–2011, and we have detailed production plans for this district as well. In this paper, correlation between temporal subsidence processes monitored using InSAR images between 2005–2013 in Bahcelievler district and coal production under this district between 2005–2011 is studied. Further, based on the obtained results, a standard vertical motion coefficient is determined for the study area.


Over time, cavities dug out during underground mining activities collapse due to the pressure of the layers above and this generates a vertical movement from the surface down to the center of the cavity (Fig. 3). This vertical movement creates the vertical deformation or subsidence on the surface.

Fig. 3. Generalized cross section of subsidence [Tandanand and Powell 1982]

Effects of vertical motion vary depending on depth of mine, size of mine, geological structure of the layers above and time. 70–80% of the subsidence is experienced on the surface during active mining and the rest is felt within 1–2 years or a little longer. The former is called active subsidence and the latter is called remaining subsidence [17].

The general formula for the maximum vertical motion that might be experienced on the surface due to mining activities underground is given by [17] as


where m is the thickness of coal extraction panel, a is the subsidence coefficient for the basin, z is time parameter and α is the slope angle of the panel. The unit of Smaks is in meters. This formula may be enhanced by adding parameter e [2] as


where e is calculated as


Based on the practical experiences gained in British mines, if the working panel width is given by W and H is mine depth and if W/H ratio is less than 1.4, it is called sub-critical production. In this case, observed vertical motion will be less than maximum vertical motion. If W/H ratio is equal to 1.4, it is called critical production and if it is greater than 1.4, it is called supercritical production/extraction. In both cases, it is expected that actual vertical motion will be close to maximum vertical motion [17] (Fig. 4). Additionally, based on the experiences gained in Indian mines, if W/H ratio is between 0.2 and 0.8, no vertical motion is experienced on the surface [15].

Fig. 4. Relationship among panel width, depth and subsidence [Sheorey et al. 2000]

The area that will be affected on the surface due to vertical motion caused by extracted seams underground can be predicted based on the limit angle of the vertical motion itself. If the coal seam is horizontal, practical experiences dictate that the limit angle will be more or less 35 degrees. Also, the maximum settlement on the surface happens perpendicular to the coal production panel. As can be seen in Figure 5, depending on the slope of the vein, limit angles and therefore the shape of the affected area on the surface will vary. Relying on these data, one can openly state that as the production level deepens, sinkage on the surface gets larger; whereas, the amount of sinkage gets lower. The radius of influence area of the vertical movement on the surface is calculated as


where R is radius of influence zone and H is depth of the production panel [17]. ξ is explained in Figure 5. The unit of R is in meters.

Fig. 5. Dip of seam and subsidence limit angles relationship [Whittekar and Reddish 1989]

In order to be able to predict the amount of vertical motion that might occur in a point inside the determined influence area, three methods; namely, empirical, quasi empirical and analytical methods have been developed. Briefly speaking, empirical and quasi-empirical methods are classified as graphical method, profile function method, influence function method, incremental design and stochastic models; whereas, analytical techniques are closed form elastic solution, numerical methods (Finite Element Method, Distinct Element Method, Finite Difference Method, Boundary Element Method) and mechanical models [10, 13].

2.1. Monitoring surface deformations by InSAR technique

SAR data is an image which consists of phase differences measured between a satellite and a ground surface. InSAR processing for deformation monitoring requires a pair of SAR data acquired at different times. One of the data pair is called the master image, and the other the slave image. The slave image is coregistered and resampled with respect to the master image. After that, comparison of the slave image to the master image yields an initial phase interferogram. Phase differences in the initial interferogram are expressed by


where Φorbit is the orbit component caused by baseline distance between two observations, Φtopo is the topographic component with respect to terrain, Φatm is a phase delay caused by reflection of microwave in water vapor layer, Φnoise is the error component caused by thermal noise, or temporal and spatial decorrelation associated with baseline distance or scattering characteristic change and Φdef represents the amount of surface deformation during the period between two observations. The components Φorbit and Φtopo are eliminated by



where Bpara  and Bperp  represent parallel and perpendicular components of the baseline, respectively. Also, h, λ, ρ and α represent elevation, wavelength, slant-range length and incidence angle, respectively [7, 12]. The component Φatm depends on total electron content of ionosphere and the parameters of the troposphere. The ionospheric part can be modeled by


where K = 40.28, fSAR is the specific SAR frequency, and the remaining term on the right-hand side that is considered as the multiplication factor is the mapping function which maps the zenith delay to the SAR incident angle Θ, using the earth radius Re and the height of the ionospheric layer hsp  [11]. The tropospheric part is calculated based on the principals of the Saastamoinen and Hopfield models [8, 14]. Φnoise is the error component caused by thermal noise, or temporal and spatial decorrelation associated with baseline distance or scattering characteristic change, and minimized by applying filtering methods. Goldstein filtering technique is the most common filtering method which is applied for InSAR processing [6]. Figure 6 illustrates the whole procedure of the Differential InSAR (DInSAR) processing summarized above.

Fig. 6. DInSAR processing procedure [Kemaldere 2011]

2.2. Monitoring Zonguldak hard coal basin

2.2.1 Previous work

As mentioned in the introduction, monitoring Zonguldak city center and its vicinity using InSAR began in 2005. In the first stage of the studies carried out until 2007, two JERS-I images taken in 1995 and sixteen RADARSAT-I images taken between 2005–2006 are utilized. For ground truthing purposes, between 2005 and 2007 parallel to the years when RADARSAT images were acquired, GPS measurements were also made. The results up to this point are published in [1].

If the first stage of the studies is summarized with reference to Bahcelievler district which is the focus area of this study, as can be seen in Figure 7, there is no motion visible in the JERS-I interferogram of 1995 in Bahcelievler district; nevertheless, there is apparent vertical motion in Kozlu, Karadon and Uzulmez under which hard coal production was taking place during this period. Here, “vertical motion” means deformation in line-of-sight; therefore, for the following sections “vertical motion” presented by interferograms should be taken as the deformation in line-of-sight.

Fig. 7. Deformed areas determined by InSAR using JERS-I images [Akcin et al. 2010]

In the analyses carried out between 2005 and 2006 using RADARSAT-1 images, a constantly evolving systematic subsidence was experienced. The amount of the subsidence growth is 5.6 cm for the period of 15 months between July 24, 2005 and October 23, 2006 (Fig. 8).

Fig. 8. Deformation obtained from RADARSAT I data for the 15 month period [Akcin et al. 2010]

2.2.2. Monitoring performed in this study

In order to be able to monitor subsidence in Bahcelievler district, until 2010 ALOS PALSAR and after 2010 TerraSAR-X data are used. The details of the interferograms are given in Table 1.

Table 1. Details of interferograms


Master Image Acquisition Date Slave Image Acquisition Date Temporal Baseline
Perpendicular Baseline
ALOS PALSAR 2006/09/24 2007/12/28 460 1923.3
ALOS PALSAR 2007/12/06 2010/03/13 828 56.4
TerraSAR-X 2011/03/15 2012/06/19 462 69.3
TerraSAR-X 2012/06/19 2012/12/12 176 41.1
TerraSAR-X 2012/12/12 2013/06/06 176 93

These images, listed in Table 1, are processed using DORIS software following the procedure outlined in Figure 6. During the processing, SRTM 90 m data is employed for topographic corrections. ALOS PALSAR satellite operates in L band and TerraSAR-X satellite operates in X band, as such, the wavelengths they use are 23.6 cm and 3.1 cm, respectively. The fringes on the interferograms which are obtained after DInSAR give the deformation on the Earth surface in the light-of-sight. Since signals travel from the satellite to the Earth and bounce back to from the Earth to the satellite, fringes depict two times the actual deformation. Thus, deformation must be divided by two. Using the said wavelengths, the fringes on the interferograms which are obtained analyzing the ALOS PALSAR data display 11.8 cm deformation in the light-of-sight and the fringes on the interferograms which are obtained analyzing the TerraSAR-X data show 1.5 cm deformation in the light-of-sight.

Interferograms rendered by processing the PALSAR images taken in September 2006 and December 2007 can be seen in Figure 9. In this 460 day period, 6 cm subsidence is determined in the area of interest, and between December 6, 2007 and March 13, 2010 subsidence had continued to grow and within the following 828 days, 18 cm subsidence is determined (Fig. 10).

Fig. 9. Deformation obtained from ALOS PALSAR data for the period between September 2006 and December 2007

Fig. 10. Deformation obtained from ALOS PALSAR data for the period between December 2007 and March 2010

Since ALOS PALSAR data were no longer available after 2010, the study was continued using X-band TerraSAR data. Therefore, six TerraSAR-X images were acquired for the dates given in Table 1. By processing these images, interferograms which are seen in Figures 11, 12 and 13 are generated.

In this respect, the interferograms below read 6.8, 3 and 1.5 cm subsidence for the periods of March 2011–June 2012 (462 days), June 2012–December 2012 (176 days) and December 2012–June 2013 (176 days), respectively.

Fig. 11. Interferogram obtained from TerraSAR-X data for the period between March 2011 and June 2012

Fig. 12. Interferogram obtained from TerraSAR-X data for the period between June 2012 and December 2012

Fig. 13. Deformation obtained from TerraSAR-X data for the period between December 2012 and June 2013


Under the area of interest where vertical motion has been experienced, five coal production panels driven at different dates are detected (Fig. 14). The production parameters and vertical movement parameters determined using these parameters and Eq. 2 and Figure 5 are listed in Table 2. Vertical motion influence areas on the surface computed based on the limit angles given in Table 2 are portrayed in Figure 14. Since driven veins are in a sense a continuation of each other, as can be seen in the figure vertical motion influence areas substantially overlap. In this regard, if a holistic solution with 50 degree average slope is pursued, the total vertical motion influence area covers all the individually determined vertical motion influence areas. Total area of the coal production panels is 54 606 m2, and calculated influence area is 1 344 165 m2. As such, using Eq. 2, amount of maximum vertical motion is estimated to be 10.9 cm.

Table 2. The production parameters and vertical movement parameters
Production Year Rise
Dip Height
Panel Area driven
Dip of seam
Seam Thickness
Coeff.  a* Limit angles
(ξ1, ξ2, ξ3)

1997–2007 -425 -560 8222 70 2.5 0.9 15, 20, 58 3.8
2003–2005 -425 -560 4860 76 2.5 0.9 15, 20, 50 3.3
2006–2007 -293 -356 3168 30 2.5 0.9 10, 25, 60 0.3
2007–2010 -300 -360 22 473 29 2.5 0.9 10, 25, 60 2.1
2010–2011 -230 -300 15 882 29 2.5 0.9 10, 25, 60 1.5
Total -230 -560 54 606 ≈50 2.5 0.9 7, 30, 65 10.9

Fig. 14. Panels driven and their influence areas on the surface (source: GoogleEarth)

Figure 15 illustrates the model of the cumulative subsidence obtained from InSAR monitoring with the driven coal panels and seams overlaid onto the satellite imagery of the interest area. Original surface of the area of interest before vertical motion is seen in Figure 16. In Figure 17, underground hard coal production panels and deformation on the surface caused by these panels is pictured. If these figures are examined carefully, it is discerned that deformations on the surface coincide exactly over the coal production panels underground. The vertical motion influence area, which is thought to be induced by the aforementioned underground coal production panels, is delineated by a white line in Figure 15. A slight surface deformation is apparent in the south-east of the area. The source of this parted deformation is thought to be other underground coal production panels or natural development of landslides.

Using Figures 15, 16 and 17, in Figure 18, south-north cross-section of subsidence progress is shown. With this figure subsidence-force relation is ascertained, and nearly all observed subsidence is within the predicted subsidence area. The area where maximum vertical motion is experienced is quite close to the estimated location. Nonetheless, compared to expected subsidence area, observed subsidence occurred in a smaller area. Based on the figure, the limit angles are 11, 28 and 31, respectively. Except dip limit angle, others are pretty close to the predicted limit angles.

In Figure 19, periodical (red line) and cumulative increase (blue line) of observed vertical motion is presented. If we look at the figure, we may deduce that observed vertical motion is accelerating between 2005–2006, steady between 2006–2008, speeding up again between 2008–2010, and afterward slowing down towards 2013. As can be seen in Figure 18, in this area coal production had continued in five panels between 1997 and 2011. If we look at InSAR analyses results, we see that during 2005–2006 deformation is increasing, during 2006–2008 no significant deformation is observed and during 2008–2010 and rapid increase in deformation is experienced and from this period onward deformation is gradually decreasing until 2013. This is the reason that we stated that if this situation is evaluated in the light of the production dates given above, it fits to the statement in literature “70–80% of the vertical motion is experienced on the surface duration of active mining and the rest is felt within 1–2 years or a little longer” [17].

Based on InSAR data analyses, the total area of the vertical movement is 353 499 m2 and maximum movement is 40.9 cm. Using observed vertical motion values and Eq. 2, if we use reverse analysis, the vertical motion coefficient “a” for the area of interest is 1.6. This value is rather high compared to the values published in the literature. Hence, in the study area there must be something else which makes deformation more severe. As the karstic formation map of the area of interest is examined, existence of intensive karstic formations is noticed beneath the area (Fig. 20). Probably, vertical motion prompted by underground mining activities triggered karstic formations close to the surface. As a consequence, observed surface deformation appears much higher than the predicted deformation.

Fig. 15. Model of the cumulative subsidence overlaid onto the satellite imagery of the interest area. The closed red line in the subset Fig. shows the estimated subsidence effect area while the white one shows the monitored subsidence (source: Google Earth)

Fig. 16. Original surface of the area of interest before vertical motion

Fig. 17. Hard coal production panels underground and deformation on the surface

Fig. 18. South-North cross-section of vertical motion evolution

Fig. 19. Periodical and cumulative increase of vertical motion

Fig. 20. Karst and fault formation of study area and its surroundings. Black dashed line indicates possible fault and blue lines display karstic caverns determined by resistivity measurements (General Directorate of Mineral Research and Exploration Zonguldak Karst and Fault Map)


InSAR monitoring of Zonguldak urban area yields maximum subsidence of 41 cm between 2005 and 2013. This result is much higher than the value of 11 cm obtained from Eq. 2. For 41 cm vertical motion, the subsidence coefficient must be 1.6 which is again much higher than the interval of 0.1–0.95 given in the literature. It is suspected that the reasons of these unusual results are shallow karstic cavities lying just below the area on the vertical motion happening on the surface.

Temporal evolution of the vertical motion displayed in Figure 19 is natural progression of the motion; that is to say, as underground coal production progressed deformations on the surface began and then as the production increased deformations became more severe, and after 2011 in the two years following the end of the production, motion has extenuated.

Observed subsidence occured within the borders of the predicted subsidence. However, because the lower end limit angle calculated using the observed vertical motion is smaller than the predicted angle of 34º, observed vertical motion area came up much smaller than estimated vertical motion area.

As a result, since in large areas it is not possible to keep all the physical parameters under control, mathematical models which are put to use for the estimation of geological movements such as subsidence, land slide etc. may not always produce results close to reality. In the example of the city of Zonguldak, the company began coal mining under the district which is largely favored and heavily inhabited in the city by stipulating that this deep coal production would cause limited deformations on the surface. Yet, existence of intensive karstic formations below the area pulled the effect of vertical motion up to dramatic levels and consequently many structures sustained damage and some of them had to be evacuated. It is solely due to the availability and accuracy of InSAR data that makes it possible to monitor these type of changes on the Earth.


  1. Akcin H., Kutoglu H.S., Kemaldere H., Deguchi T., Koksal E., 2010. Monitoring subsidence effects in the urban area of Zonguldak Hardcoal Basin of Turkey by InSAR-GIS integration, Nat. Haz. & Earth Sys. Sci., 10 (9), 1807–1814.
  2. Bischoff W., Bramann H., Dürrer F., 2010. Das kleine Bergbaulexikon, 9. Edition, Vgë Verlag GmbH, 393 p.
  3. Can E., Kuşcu T.,Mekik C., 2012. Determination of underground mining induced displacements using GPS observations in Zonguldak-Kozlu Hard Coal Basin, International Journal of Coal Geology, 89 (1), 62–69.
  4. Deguchi T., Kato M., Akcin H., Kutoglu H.S., 2006. Automatic Processing of Interferometric SAR and accuracy of surface deformation measurement, Proceedings of SPIE – The International Society for Optical Engineering, Sweden, 6363, 636309.
  5. Deguchi T., Kato M., Akcin H., Kutoglu H.S., 2007. Monitoring of mining induced land subsidence using L- and C-band SAR interferometry, International Geoscience and Remote Sensing Symposium (IGARSS), 4423253, 2122–2125.
  6. Goldstein R.M., Werner C.L., 1998. Radar interferogram filtering for geophysical applications, Geophysical Research Letters, 25(21), 4035–4038.
  7. Hanssen R.F., 2001. Radar interferometry – data interpretation and error analysis, Springer, New York.
  8. Hopfield H., 1971. Tropospheric Range Error at Zenith, Committee on Space Research, 14th Plenary Meeting, working group 1, id. number a.15, edited by Applied Physics Laboratory, The Johns Hopkins University, Maryland.
  9. Kemaldere H., 2011. Monitoring of undercity mining and subsidence effect by differential InSAR technique in Zonguldak Metropolitan Area, Ph.D. thesis, Zonguldak Karaelmas University (In Turkish).
  10. Kratzsch H., 1983. Mining Subsidence Engineering, Springer – Verlag, 41, 153.
  11. Likoka S.A., Karathanassi V., 2007. A statistical analysis of the error produced by the tropospheric path delay calculation, Proceeding of Envisat Symposium 2007, 23–27 April 2007, Montreux, Switzerland.
  12. Massonnet D., Feigl K.L., 1998. Radar interferometry and its applications to changes in the earth’s surface, Rev. Geophys., 36(4), 441–500.
  13. Peng S.S., 1992. Surface subsidence engineering: Littleton, Colorado, Society for Mining, Metallurgy, and Exploration, Inc., 161 p.
  14. Saastamoinen J., 1973. Contribution to the Theory of Atmospheric Refraction, Bulletin Geodesique, 105–107.
  15. Sheorey P.R., Loui J.P., Singh K.B., Singh S.K., 2000. Ground subsidence observations and a modified influence function method for complete subsidence prediction, Int. J. Rock Mech. Min. Sci., 37, 801–818.
  16. Tandanand S., Powell L.R., 1982. Assessment of subsidence data from the northern appalachian basin for subsidence prediction, Report of Investigations – United States, Bureau of Mines.
  17. Whittaker B.N., Reddish D.J., 1989. Subsidence: Occurence, Prediction And Control, Elsevier, Amsterdam, The Netherlands, 528 p.

Received: 18.09.2018
Reviewed: 13.03.2019
Accepted: 24.04.2019

H.S. Kutoglu
Geomatics Engineering Department, Faculty of Engineering, Bulent Ecevit University, Zonguldak, Turkey

email: kutogluh@hotmail.com

M. Berber
Department of Civil and Geomatics Engineering, California State University, Fresno, CA, USA

email: muberber@csufresno.edu

H. Kemaldere
Geomatics Engineering Department, Faculty of Engineering, Bulent Ecevit University, Zonguldak, Turkey

email: kemaldere@hotmail.com

Responses to this article, comments are invited and should be submitted within three months of the publication of the article. If accepted for publication, they will be published in the chapter headed 'Discussions' and hyperlinked to the article.