Assessment of geosites as a basis for geotourism development: A case study of the Toplica District, Serbia

2.1 Branko’s tower (GS-1)

The spectacular outcrop of columnar andesite basalt is located about 12 km south of Kuršumlija, about 1.5 km southwest of the village of Rudare, at 747 m a.s.l. (Figure 2(1)). It is located in the protected zone of the Radan Nature Park. The andesite-basalt pillars resemble a wall about 30 m long and up to 3 m high. This formation, which runs roughly in the direction of NNE-SSW, consists mainly of five-sided pillars, usually about 30 cm in diameter. The pillars are composed of hornblende-pyroxene andesite. The product is the second volcanic phase of the large Lece Magmatic Complex (LMC). The structure is actually a shallow dike-like intrusion. The columns are separated by joints or fractures that generally form when the lava shrinks, usually during cooling. They cut through the volcanic mound of hornblende andesite volcaniclastics, which are products of the first volcanic phase [33].

Figure 2

(1) Spectacular andesite-basalt stone pillars at Branko’s tower; (2) Flysh rock formation at the village of Krčmare; (3) Limestone massif and cave Treska; (4) Paleovolcanic dome and archeological site Ivan’s tower; (5) Cone-shaped formation of volcanic stone pillars Devil’s City; (6) Middle Miocene paleontological site Prebreza.

2.2 Landform group in flysch rocks near the village of Krčmare (GS-2)

Flysch is a dense marine sequence of turbidites and gravity mass flows embedded in the overthrust and fold belt of a growing orogenic belt [34]. To be recognized as flysch, the deposits must have appropriate lithologic and depositional characteristics. In the first case, it is the gradual deposition of suspended sediments from deep-water turbidites to shallow-water shales (clays and silts) and sandstones. At the base of the flysch sequence are coarse-grained sediments, breccias, and conglomerates. Differential erosion of flysch sediments, especially in the strata deviating from the main topographic surfaces, leads to the formation of very complex landscapes with specific landforms such as rocks in the form of pillars (stacks), cones, tables, mushrooms, haystacks, abri-scarps escarpments and rock shelters, isolated ridges, and boulder fields [35]. Vivid examples of reliefs created by selective weathering of Upper Cretaceous flysch breccias and conglomerates within the so-called Toplica flysch [30] are recognized in the village of Krčmare (right bank of the Toplica river, about 4 km from Kuršumlija). The figures are locally known as “Mara’s petrified stacks.” The most vivid are Šiljak (sharp stone, a spike), Čučavac (squatting), and Kamalj (the stone) (Figure 2(2)), but there are also a number of head-like rock formations up to 1.5 m in diameter. It is carved into well-cemented conglomerates and is about 6 m high. Čučavac is a 3.5 m high rock without a base, which looks like a Dutch clog. The third relief figure, Kamalj, resembles a lying fossilized prehistoric animal.

2.3 Treska limestone formation (GS-3)

Treska limestone formation is located about 50 km northwest of Kuršumlija and about 20 km west of Lukovska Spa, between the villages of Seoce, Štava, and Trešnjica. Treska (1.439 m a.s.l.) consists of tectonized massive and sandy Triassic limestones (Figure 2(3)). It is one of the dominant peaks on the eastern slopes of the Kopaonik Mountains. Its appearance in the form of a double tooth makes it stand out from other mountain peaks. It is accessible only from the west side, and the ascent from the village of Štava is about 7 km long. It is characterized by difficult terrain, but also by exceptionally preserved biodiversity. It is especially attractive for mountaineers and lovers of geotourism.

The almost unknown and little explored Treska Cave is located at the foot of the southeastern part of the limestone massif. The estimated explored area of the cave is about 0.23 km², which is about 40% of its total area. It is relatively poor in speleothems. Preliminary investigations of the cave show that it consists of 11 galleries, but they are narrow and inaccessible. The entrance to the cave is now collapsed, so the cave remains unexplored and its total length is unknown.

2.4 Ivan’s tower (GS-4)

Ivan’s tower is a paleovolcanic dome (886 m a.s.l.) and a partially explored archeological site (Figure 2(4)). It is located on the western slopes of the Radan Mountains, 33 km south of Kuršumlija, in the immediate vicinity of the villages of Zagrađe and Ivan Kula. The paleovolcanic caldera is part of the LMC. The preserved archeological remains are related to a medieval fortress built on top of a paleovolcanic dome with a rectangular base and a height of 14 m. The area of 65.47 ha around Ivan’s tower is protected by the second level of protection within the special geo(morpho)logical nature reserve “Kosanica-Đavolja varoš” [36]. The archeological remains of the medieval fortress have been protected by decree as a cultural monument since 1956 [37].

2.5 The Devil’s City geocomplex (GS-5)

The Devil’s City geocomplex is the best-preserved part of the LMC. The area of about 4.3 km² has been under protection since 1959. In 1995, it was declared as a Natural Monument of protection category I, and since 2009, it has been on the Tentative list of UNESCO. Besides the unique geomorphological phenomenon of earth figures, it also includes several different geosites (Figure 2(5)). It is located at an altitude of 700–720 m a.s.l., on Radan Mountains, about 29 km southeast of Kuršumlija.

The Devil’s City is the most important geomorphological natural monument in Serbia, formed by erosion [38]. Striking erosional forms (a total of 202 earth/stone pillars, earth columns, soil/clay pyramids) are located between two gorges, the Devil’s Gorge and the Hell’s Gorge. Pyroclastic rocks are mostly represented by weakly bonded volcanic breccias and tuffs. These rocks are exposed as horizons ranging in thickness from a few meters to 100 m. Notable features of the Devil’s City are the predominance of pyroclastic rocks over andesites and intense hydrothermal activity. Cone-shaped forms of poorly sorted volcanic material with large andesite caps reach heights of 2–15 m, widths of 4–6 m at the base, and heights of 1–2 m at the top. These ephemeral forms disintegrate relatively quickly when they lose their protective “cap” of andesite blocks [39].

Devil’s City is also characterized by rare hydrological phenomena, springs with rare hydrochemical properties related to the last phases of volcanism. The first spring known as “Devil’s Water” is extremely cold, acidic, and highly mineralized water [40] with radioactive content [41]. The second spring “Crveno Vrelo (Red Spring)” is less acidic but contains heavy metals, especially iron. The wider area of Devil’s City is rich in ores (iron, aluminum, gold, and silver), and was part of a well-known mining area during the Middle Ages. The presence of several mining shafts are archaeological remains of mining from the thirteenth century.

2.6 Paleontological site of Prebreza (GS-6)

Several fossil sites exist in the Toplica district. The paleontological site “Prebreza” is one of the best preserved and most important European sites of vertebrates from the Middle Miocene (Figure 2(6)). It is located on the southwestern slopes of the Veliki Jastrebac Mountain, in the northwestern part of the Tertiary basin of Toplica, near Blace. The area has been protected since 1963. In 2009, it was declared as a special nature reserve and is protected by the Inventory of Geological Heritage of Serbia. The fossil remains of Prebreza are dated to the middle Miocene Badenian–Sarmatian [42]. The site was discovered by erosion of the coastal slopes of the Gluvi Potok (right tributary of the Blatašnica River). It extends on both sides of the watercourse in the form of a subvertical profile several hundred meters long. The geological structure of the Prebreza area consists of sandstones and clay sandstones with mica, covering pelite and tuff with coal.

3.1 Data collection

In the initial phase, available literature and spatial planes were used as a source of spatial data. Digitized and georeferenced geologic maps 1:100,000 [30]; geomorphologic maps 1:300,000 [43], and military topographic maps 1:25,000 and 1:50,000 were used. By selecting different sources, we searched for symbols representing geosite features (cliffs, fractures, gorges). Geographic coordinates are defined by determined geosites. Available data via Google Earth and OpenStreetMaps applications in the research area were searched to create the preliminary GIS database. The focus was on finding photos of specific natural formations, locations, and viewpoints that have already been marked by climbers and geotourists. The free online database of the Information system of the immovable cultural property was also used (https://nasledje.gov.rs/index.cfm?jezik=Serbian_CIR) [44]. The location of geosites was marked in the field using GPS and mobile applications.

To reduce the number of sites to be checked in the field, pre-filtering was done based on geological significance, attractiveness, and size. Some of the geosites found proved to be less important, difficult to access or inaccessible, and of small size. Such sites were excluded from the GAM. In addition, the number was reduced by grouping them by location or similarity of origin. Finally, due to the complexity of the survey, six geosites were selected for analysis. A detailed report of the survey of geologic features was written for each geosite. Photographs of each geosite visited were taken.

3.2 Questionnaire designs

The application of the GAM was preceded by a survey in which a total of 52 people participated. The geospatial experts were carefully selected for their expertise in geography, geology, geoheritage, and geotourism. They are familiar with the research topic and the study area. The anonymity of the questionnaire preserved the credibility of the responses. The research was conducted in February and March 2023. Statistical calculation and visualization of the results were performed using SPSS 20 (Statistical Software for Social Sciences) and MC Excel. In preparing the questionnaire, the purpose of the research was first explained and each question and group of answers was simplified individually. The questionnaire is divided into two parts. The first part contains questions about the socio-demographic profiles of the participants (gender, age, place of residence, and education), as shown in Table 1. An important question for the study was related to the type of transport most often used on the journey to and during the stay in destinations where the participants could choose between different modes of transportation. In the second part of the questionnaire, each participant rated each of the 27 sub-indicators. A five-point Likert scale ranging from 0.00 to 1.00 was used for the rating (with 0.00 representing unimportant and 1.00 representing very important). The survey was conducted online, and the Google Forms application was used to design the questionnaire. The link to the questionnaire was sent to each participant’s email address and can be found in the supplementary material section of the study. The aim was to have only geo-experts participate in the survey, thus ensuring the validity of the sample.

The questions included in the survey have already been tabulated in a large number of studies that use the GAM and/or M-GAM methods. Therefore, they are not tabulated again in this study. The values obtained by statistical calculations based on surveys and the final results of the geosite assessment are presented.

3.3 Geosite assessment model (GAM)

In order to present and quantitatively evaluate the macro and micro forms of relief in the Toplica district, the GAM developed by Vujičić et al. [24] was chosen. It is one of the most complex and well-defined methods. Numerical evaluation is based on surveys and field research. The modified version (M-GAM), developed by Tomić and Božić [25], is preceded by identical survey and field research. This method includes the evaluation and opinion of tourists, which often and decisively affects the final results of the assessment. Both models, through numerous studies, have already been applied in the evaluation of different areas of Serbia and neighboring countries [24,25,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60]. Most of the research geosites of the Toplica district have not been publicly exposed, have not yet been geotouristically evaluated, and rarely been visited individually, so the GAM method proved to be adequate for this research.

In the analysis using the method GAM, two groups of values are evaluated. The first group, Main Scientific Values (MV), was created according to the natural characteristics of geosites. They consist of three groups of indicators – Scientific/educational value (VSE), Scenic/aesthetic value (VSA), and Protection value (VPr), with a total of 12 sub-indicators evenly distributed. For each geosite and each of the sub-indicators, the geospatial experts score one of the following numerical values: 0.00, 0.25, 0.50, 0.75, or 1.00. The total summary of these three groups of indicators is the final score for the MV (equation (1)).

where SIMV _i represents the sub-indicators of the MV (i = 1, …, 12).

The second group, the Additional Values (AV), refers to the evaluation of the degree of a tourist-functional position of the geosites. They consist of two groups of indicators – Functional values (VFn) and Tourist values (VTr) – with a total of 15 sub-indicators, unevenly distributed into six and nine sub-indicators. Following the same principle, the same summation procedure is performed to obtain AV (equation (2)).

where SIAV _i represent the sub-indicators of the AV (i = 1, …, 15).

According to the definition of the method GAM, the addition of the two groups of values gives the total result for each geosite (equation (3)).

A matrix divided into a total of nine fields was created for the graphical representation of the overall evaluation results. The H axis consists of the MV score, ranging from 1 to 12. The Y axis consists of the AV score, ranging from 1 to 15. Depending on the final result, a specific field of the matrix and the so-called “tourist value” is determined for each geosite. In order to better understand the values obtained, the standard classification (low, moderate, high) has been slightly modified. The level “moderate” is divided into “moderately high” and “moderately low” to highlight the tendency of geosites to have a high or low tourist value.

3.4 Data analysis

The remaining analyses were performed using the GIS. GIS and numerical modeling are becoming powerful tools [61] for inputting, editing, analyzing, creating, and improving spatial data. GIS Software package Quantum Geographical Information System 3.18. Zurich (QGIS) was used for data processing. For spatial analysis, the System for Automated Geoscientific Analyzes extension was used in the software. The methods of digitization, vectorization, and geo-referencing are used [27]. A raster grid with geo-referenced coordinates is downloaded. Then, the boundaries of the Toplica district were determined in vector format and connected to the locations of geosites within the district. The SRTM Digital Elevation Model (DEM) was also used. Using the GIS software, the boundaries of the Toplica district were cropped for processing purposes. Elevation zones were determined using QGIS within DEM. The IDW interpolation model was also used for geostatistical analysis and data modeling. Surface grids were created based on the defined geosite locations as well as the data obtained through the analysis of GAM. Using QGIS, the areas of the gravitational influence of the studied geosites were measured and analyzed based on the created raster representation.

Although there are several methods, the IDW method is preferred because it assumes that each measured point influences the locations that decrease with distance. This is one of the most commonly used deterministic interpolation methods for creating analytical distribution maps. It is considered a suitable method for the distribution of grouped points (clusters) in local space, such as the geosites under study. IDW is one of the spatial interpolation techniques that allows conversion of discrete measurements into a continuous spatial distribution. It is based on the assumption that neighboring values contribute more to the interpolated values than distant observations. For predicting values at any unmeasured location, IDW uses the measured values in the vicinity of the prediction location. It is based on the assumption that neighboring values contribute more to the interpolated values than distant observations. For predicting values at any unmeasured location, IDW uses the measured values in the vicinity of the prediction location. It is characterized by ease of use, clear and efficient interpretation of results [62], intuitiveness, and the “bull’s eye” effect. The IDW interpolation method has not been used in the analysis of GAM. In the field of tourism, the method has been used with the implementation of GIS in the interpretation of the GIS-based tourist index [63], the analysis of social networks for tourists [64], the comparison of tourist images of the destination [65], the evaluation of the tourism potential of rural areas [66], the evaluation of the tourism suitability of wetlands [67], emotion mapping in urban areas [68], and sustainable tourism planning [69].

The general formula for IDW is as follows [70]:

where Z represents the interpolated value; z i is a known value; n is the total number of known values used in the interpolation; d i is the distance between known and interpolated values; p is the power parameter where the weighting coefficient decreases as the distance between the interpolated points increases; the power parameter as an exponent of the distance between points, is defined by the most commonly used value (2) [71]. During the research, the calculated MV and AV of each geosite are used to convert into spatial values based on the IDW interpolation method. The overall results of GAM give a better overview of the tourist value of the place.

In the second part of the study, network analysis was used to analyze the traffic accessibility of geosites by determining their temporal range. Isochronous maps (“time-contour maps” or “travel-time maps”) are based on the assumption that accessibility is a criterion that decisively influences tourist suitability. The criterion that primarily influences the tourist suitability of a geosite is its accessibility. The management of geotourism resources and the tourist development of geotopes are related to this criterion. Inaccessible sites require much more investment in infrastructure, which affects their tourism profitability. Due to the temporal distance, isochronous maps show the availability and analyze the space and points of interest (tourist infrastructure, tourist attractions, and geosites). The isochronous method of GIS has been applied in solving tourist traffic problems, determining traffic accessibility of attractions using GIS in Ukraine [72], analyzing the results of traffic accessibility of public transport in Warsaw, Poland [73], understanding the potential for the development of transit-oriented tourism [74], and analyzing the availability of roads to tourist destinations [75]. This method is also used in the multidimensional representation of natural and cultural heritage in Slovenia [76], visualization of tourism planning patterns [77], and identification of the most suitable mine sites for geotourism planning [78].

In the creation of the isochron maps, the geospatial data of the geosites in vector format and the road network of the OpenStreetMaps platform in vector format were used in one of the ORC (OpenRouteService) algorithms of the QGIS plugin.

Through the algorithm of the installed plugin, a set of polygons was created based on the road network around the given location. The measured polygons represent areas that can be reached within a certain distance in time. Travel times are calculated for each segment using average speed limits for different road types and considering different gradients or road surfaces. The parameters of the analysis were based on the results of the survey, which showed that more than 83% of the respondents use a car for tourist visits. The analyzed space has a local character; time intervals of 15 min were used. Geosites are evaluated by availability, from highly available to low available and unavailable (or slightly available), with moderately-high and moderately-low transition variants. The objective is to use a vector traffic network over a time radius to determine the traffic accessibility of geosites and analyze their suitability for tourist visits. In this way, it will be shown how isochron maps can be a gateway for more efficient geotourism organization and planning in the Toplica district.

4.1 GAM numerical analysis

The previously described geosites of the Toplica district were evaluated using the GAM methodology. The results (Table 3) show the uniformity of indicator scores for the MV and AV of all geosites. The results also indicate a large mutual difference in the total GAM values between them. Such a relationship indicates a great potential but uneven tourism development of the studied area.

The ratings of the MV indicators presented in Table 3 highlight the Devil’s City geosite. Undoubtedly, the best evaluated, with a high overall score (9.12), stands out with extremely high scores for Scientific/educational value (3.27). The consequence of sustainable tourism management is high scores for the Scenic/aesthetics and Protection sub-indicators (Table 2). Devil’s City has a relatively small area, and the structure of the relief is limited, so it is not surprising that it is rated lowest in terms of Viewpoints.

Approximately equal ratings of all Main indicators, with better rated Scenic/aesthetic sub-indicators, are the result of Prebreza’s geosite assessment. The limitation of the relief structure in this case also affects the low scores of the sub-indicators for Viewpoint. On the other hand, the highly rated sub-indicator Environmental fitting of the sites is a result of protective measures. Three geosites, Ivan’s tower (5.36), Treska (4.81), and Branko’s tower (4.32), were given lower MV scores. These sites are rated better in terms of the Scenic/aesthetic sub-indicators. In relation to the Protection indicator, they were rated significantly lower. Thus, these are sites of exceptional aesthetic and landscape value, but they are poorly known to the public and only partially studied scientifically. This is confirmed by the extremely high rating of the Environmental fitting of Ivan’s tower geosite. The lowest rated geosite is GS-2, with a low overall rating of the MV (3.65). The very low scores of the indicator of Environmental fitting of the sites and especially the extremely low scores related to Protection are evidence of its unrecognizability.

The AV contain important indicators of geotourism development. As can be seen from Table 3, Devil’s City also stands out here with a high rating (9.76). This rating is a consequence of the high-rated indicators of Tourist values and VFn, especially the sub-indicator of Promotion. The geosite of Prebreza is in second place with an overall score of 6.57. The proximity of the municipal center of Blace has a significant impact on the higher rating of the Tourist values (3.86) compared to the VFn (2.71). This relationship also applies to the sub-indicators. The Accessibility and Road network sub-indicators were rated the highest, while the Number of visitors and Guide service were rated the lowest. The ratings indicate a high geotourism potential, but also its low utilization. The other geosites GS-1 to GS-4 are rated low in terms of AV. The ratings of all geosites are between 3 and 4, indicating low tourism recognizability and modest infrastructural investments.

After comparing the final results, the differences between the MV and AV become clearer through the positions of the 6 geosites in the GAM matrix (Figure 3). The best-valued geosite GS-5 Devil’s City is also the best positioned, in the high-ranking field Z32 of the GAM matrix. The position indicates highly evaluated MV and moderately high evaluated AV. In total, the tourist value of this geosite is 18.88 (Table 3). In second place, according to the total score of GAM, is the site GS-6 Prebreza. The MV of moderately high level and the AV of moderately low level, position this geosite in the field Z22 of the GAM matrix. The overall tourism value, defined as moderately low, is scored as 13.21. The three geosites, GS-1, GS-3, and GS-4 are positioned in field Z21 of the matrix. Their MV are defined as moderately low, and AV are defined as low, which significantly affects their position in the GAM matrix. The overall rating of all geosites is about the same. The scores of the indicator groups of the geosite GS-4 Ivan’s tower are slightly higher than the scores of the other geosites. Thus, it is best positioned in the Z21 area, the GAM matrix. The worst-rated geosite GS-2 is located in the Z11 field of the GAM matrix and has the lowest ratings of the indicator groups. Its tourist value can be defined as low.

Figure 3

Scatterplot diagram of the GAM results – Branko’s tower (GS-1); Landform group near the village of Krčmare (GS-2); Treska limestone formation (GS-3); Ivan’s tower (GS-4); Devil’s City (GS-5); Prebreza paleontological site (GS-6).

4.2 GAM geostatistical analysis

Based on the results of the numerical analysis of GAM, the values were interpolated into the surface raster. By visualizing the GAM indicator groups, a comparative distribution map was created (Figure 4). Since QGIS software also provides the ability to measure surfaces, it was used to analyze the gravitational influence of geosites. The results of measuring the gravitational influence of geosites for MV, AV, and final GAM Tourist Values are presented in Table 4. Comparative examination of the interpolated areas in Figure 4 shows that both indicator groups almost equally participated in most of the geosites. There is a balance between scientific-educational and tourist-functional indicators.

Figure 4

Distribution map of GAM values: Left – GAM MV; Right – GAM AV.

According to Figure 4, two gravity centers are defined. They are located on different (opposite) parts of the interpolated area. It is more meaningful in terms of geosite Devil’s City. The gravity range of the MV for this geosite is more expressive but has a smaller range of influence. As a result, there is a greater deviation in its ratings compared to the ratings of the other geosites. The measurement showed that the area is 224 km² (Table 4). The gravitational area for the AV of this geosite is much larger (520 km²), but less meaningful because it has a smaller deviation of the GAM ratings. Another milder and less meaningful centroid is associated with the area of the gravitational influence of the geosite Prebreza. As a result of the smaller deviation of its ratings concerning the ratings of the neighboring geosites, its gravitational area for MV is smaller (171 km²) – in terms of values, the value of GAM is on the lower limit of the moderately high level). The area of influence for AV is larger (364 km²) because the scores of neighboring sites for this group of indicators are very low.

The deviation of the interpolated areas can be seen on Ivan’s tower, as the values of the Main indicators are significantly higher than the values of the Additional indicators. The gravity area of influence is larger for the MV (95 km²) than for the AV (40 km²). This relationship is a consequence of the large differences in the GAM values (especially the AV) of these two neighboring geosites. The areas of the gravitational influence of other geosites overlap, so it can be assumed that they jointly participate in the interpolated area. This is particularly clear for the AV, where the three geosites GS-1, GS-2, and GS-3 form a common area between the two gravitational centers due to their low-scoring indicators. For the MV, the geosites Branko’s tower and Treska, defined by the same distance segment, form a common gravitational influence area of 1,307 km². The worst evaluated gravity influence area of the geosite GS-2 is a single area of 108 km².

Based on Table 3, the total values of GAM were interpolated into the distribution map (Figure 5). On the map, separate areas of the geosites GS-1, GS-2, and GS-3 can be seen, then the common gravity area of the geosites GS-4 and GS-6, and the strongly pronounced gravity center of the geosite Devil’s City. This ratio results from the total values of GAM (Table 3) and is a more plastic representation of the ratio in the GAM matrix. As a consequence of the best positioning in the GAM matrix, the area of gravitational influence of the geosite Devil's City (112 km²) is related to the highest segments of the range (Figure 5). The common area of geosites GS-1 and GS-2 (337 km²), and GS-3 (172 km²), are the result of the lowest segment of the range – low defined tourist value (<9.00). Geosite GS-4 and GS-6, which are better positioned in the matrix GAM, belong to the moderately low segment of the range and form a common gravitational area (1,426 km²).

Figure 5

Distribution of the final GAM tourist values.

Figure 6

Perception of the tourist accessibility of geosites in the Toplica District.

4.3 Isochron analysis

Since most of the geosites gravitate to the area of the municipality center of Kuršumlija, it can be considered the carrier of geotourism development in the district. It was chosen as the starting point for the second part of the research. As a prediction for the improvement of the tourist offer of the study area, the isochron analytical map evaluates the geotourism potential of the sites through accessibility as the initial factor of geotourism activation. Similar analytical maps could be created for other municipal centers and tourist resorts of the Toplica district in order to expand the tourist offer and create geotourism routes.

The isochrone map based on the GIS algorithm was created in the area (Figure 6). It shows that the area of the Toplica district was divided into five differently colored isochron zones. The transportation and tourist accessibility of each of the six geosites is defined by its position in one of the isochron zones. Only one geosite was located within the first isochron zone. In the immediate vicinity of Kuršumlija, the Krčmare landform group is well below the maximum time interval of this zone (maximum 15 min). This site makes it highly accessible for tourists. On the border between the first and second isochrone zones, but with the predominant belonging to the second zone, there is also a geosite, a Branko’s tower. The location accessibility, determined by a maximum radius of 30 min, is defined as moderately high. The third isochrone zone includes two geosites, Devil’s City and Prebreza. Although they are located in opposite parts of the district, their distance from Kuršumlija is almost the same. The temporal distance of not more than 45 min makes them moderately low accessible. According to the isochrone analysis, the most distant sites Treska and Ivan’s tower are included in the fourth isochrone zone with a maximum radius of 60 min. Equally distant from Kuršumlija and located at different ends of the Toplica district, their location could be classified as peripheral and their accessibility as low.