Nevşehir Castle Region in Turkey Interpreted by the Use of Seismic Surface Wave and Electrical Resistance Measurements Together

Abstract The underground city beneath the Nevşehir Castle located in the middle of Cappadocia region in Turkey with approximately cone shape is investigated by jointly utilizing the modern geophysical techniques of seismic surface waves and electrical resistivity. The systematic void structure under the Nevşehir Castle of Cappadocia, which is known to have widespread underground cities, is studied by the use of 33 separate two-dimensional profiles ~4-km long where electrical resistivities and seismic surface waves are concurrently measured. Seismic surface wave measurements are inverted to establish the shear-wave velocity distribution while resistivity measurements are inverted to resolve the resistivity distribution. Several high-resistivity anomalies with a depth range 8-20 m point to a systematic void structure beneath the Nevşehir Castle. We were able to effectively isolate the void structure from the embedding structure since the currently employed resistivity instrument has provided us high resolution quality measurements. Associated with the high resistivity anomalies there exist low-velocity depth zones acquired from the surface wave inversions also pointing to a systematic void structure where three-dimensional visualization techniques are used to show the extension of the void structure under the studied area.


Underground Cities in Cappadocia
Erciyes, Hasandağ and Melendiz volcanic eruptions dating back to Upper Miocene to Quaternary times resulted tuff layers which were then eroded by wind and water regime creating remarkable landscape famously known as Fairy Chimneys in central Turkey [23,24]. These volcanic rock layers extensive in the region were carved by early settlers to built the underground cities ( Fig. 2 -upper row). Researchers have reported that the rock settlements in Cappadocia were either totally located beneath the surface or semiunderground or carved into a cliff [25]. In Derinkuyu underground city, there are eight floors going more than 80 meters deep down into the earth. This underground city with a number of rock rooms, religious centers, store rooms, winery, kitchens and stables for live stocks is estimated to have enough space for 15,000 -20,000 people ( Fig. 2 -lower row). There are more than a hundred underground settlements in Cappadocia. Alongside the map of the Anatolian plate Fig. 3 shows locations of main underground cities discovered in Nevşehir province of Cappadocia [26]. The Nevşehir Castle region indicated by a green color circle in Fig. 3 is located close to the Nevşehir city center. The newly discovered underground settlement reported herein appears to be the center of all other underground settlements known in the region. Our field team and also people from Nevşehir city administration have collected field data from numerous locations of the underground city studied herein. A circular stone cut-and-shaped in place used to prevent the entrance into the secretive room provided by a passage big enough for only one person (upper left), another circular stone attached to a crop grinding machine (upper right), an underground room carved with supports provided by vertical columns (lower left) and another underground room having outside access through a window (lower right) are shown for some of the findings (Fig. 4), which are actully structures commonly found in other underground cities in Cappadocia.

Brief Geology
The study area is located in Cappadocia which is generally underlain by volcanic rocks of Neogene-Quaternary period. These volcanic units are generally of dirty white, grey and pink colors and of pyroclastic character as formed by the neotectonic convergent volcanism between the Eurasian and Afro-Arabian plates occurring in the eastern Mediterranean [25,26].
A volcano-sedimentary layer (Miocene-Pliocene), which is called the Ürgüp formation comprising 10 ignimbrite units (Kavak, Zelve, Sarimaden, Cemilköy, Tahar, Gördeles, Sofular, Kızılkaya, Valibaba and Kumtepe), covers most of the study area outcropping with considerable thickness while hosting interesting surface features after water and wind erosions (see Fig. 5). The Ürgüp formation with a volume of at least 1000 km3 and an area of about 40,000 km2 lithologically shows variations in both vertical and horizontal directions. Among this variation the Kavak unit (10-150 m thickness) exposed in Kavak village, which is made of intercalation between volcanic ash-flow deposits and volcanoclastics, constitutes the lowest level layer in the sequence [27]. All the ignimbrite units taking their names from the nearby villages in the Ürgüp formation are well discussed in [24] and the references therein.

Materials and Methods:
The field study lasted approximately six months. The resistivity measurements were taken in the last months of 2012 and the seismic measurements were taken in the first months of 2013. During the resistivity measurements we waited for dry weather to avoid rain and snow fall. Fig. 6 shows the profiles (~4km long -blue lines) where the profile names are in Turkish taken from street names along which the profiles are measured. In order to avoid curvatures along our profiles and also field obstacles, some profiles were measured in relatively smaller pieces. The streets on which the geophysical measurements were taken were covered by 10-cm thick concrete bricks below which ~50-cm thick sand and gravel layer took place. Under these circumtances we tried to collect the data, but the quality was low. Therefore, the brick and the sand-gravel layer on all measurement profiles were scraped off by an excavator creating ~70-cm deep trench after which we were able to collect quality data placing the sensors on the trench bottom.
The study area is close to the city center. The seismic noise level due to the cultural activities was relatively high. In order to overcome this noise condition, we designed a powerful accelerated weight drop system inspired from the system given in [28]. The machine (see middle panel in Fig. 7) was built by Nevşehir Technical and Industrial Vocational School. It was operated by two people and had 60-kg mass easily lifted by a hand-controlled system 50-cm high above the ground against two springs each 90-cm long and then released for a strong hit on the ground. In each seismic experiment we stacked 8 to 15 hits to increase the signal quality and the accelerated weight drop system was used around 1500 times including the test profiles.
A multi-channel and multi-electrode system with 5-m electrode spacing was used to conduct the ERI survey. The ERI equipment is ZZ Geo FlashRes320 resistivity instrument [29], which is a 317-channel 128-electrode, full-wave-shaped electrical data collection system and is free from any configuration. Herein we do not repeat the information regarding the ERI instrument and also the inversion methodolgy of the resistivity data since these are extensively reported in [29]. The ERI system (currently 32-electrode) does not use a traditional data collection technique (e.g., Wenner, Schlumberger, or dipole-dipole arrays) and the maximum combination of ABMN (current and potential electrodes) permits to collect a large amount of data simultaneously with 317channel system. Therefore, a large amount of data collected significantly improved the resolution of the electrical resistivity inversion.
A 24-geophone seismograph system (DAQlink III, see left panel in Fig. 7) with 4.5-Hz vertical-component geophones was used to collect the seismic data (see right panel in Fig. 7 for a sample recording). In general, the geophone spacing was set to 2-m. In some cases, depending on the field conditions the geophone spacing was reduced to 1.5-m. The shot point was placed at an offset of 30-m from the first geophone and again depending on the field conditions the source offset was occasionally reduced down to 25-m. In order to invert the one-dimensional (1-D) shear-wave velocities beneath a geophone spread we utilized the MASW technique first introduced by researchers as an extention of spectral analysis of surface waves (SASW) for 1-D shear-wave velocity structure [30][31][32]. The method directly converts timedistance domain waveform data into an image of phase velocity versus frequency and utilizes the dispersion property of surface waves for the shear-wave velocity profiling in 1-D (depth) or two-dimension (2-D, depth and shifted surface location). The data processing of the MASW method generally consists of three stages: 1) measurement of seismic surface waves generated by the seismic source, 2) extraction of the fundamental-mode phase velocity dispersion curve and 3) inversion of the dispersion curve to obtain the 1-D velocity-depth profile [30]. Both source and receivers move up along the survey line to interpret the velocity structure in 2-D. We currently moved the source and receiver array by one receiver-spread length without an overlap.

Surface Wave Analysis
The "Computer Programs in Seismology" by a researcher was utilized for the second and third MASW steps above [33]. We only inverted the fundamentalmode dispersion curve for the velocity structure. The three test models shown in Fig. 8 (left panel) provided us some insight into possible phase velocity curves (right panel) expected for the shallow structure beneath the Nevşehir Castle region. The velocity structure depicted by the black line represents a regular velocity structure steadily increasing down to the half-space at ~20-m depth. The second velocity structure depicted by the blue line has a low velocity zone at ~13-m depth and is the slowest model among the three models. The third model depicted by the red line also has a low velocity zone, but deeper at ~23-m depth. The second and third models has a deeper half-space around 33-m depth. The low velocity zones in these two models are representative of voids at different depths. The three set of phase velocity curves in the period range 0.03-0.3 s in the right panel present distinct properties to discriminate between the three velocity structures in terms of the fundamental mode considered herein. It is noted that the fundamental-mode dispersion curve corresponding to the model depicted by the red color line tends to flatten around 0.06-s period because of the low velocity zone in the model. This is also true for the blue color model, but starts at a relatively shorter period. Fig. 8 suggests that the fundamental mode phase velocity curve is good enough to delineate the low velocity zone due to a void structure at some depth.

Figure 8:
Three test models on the right and the corresponding Rayleigh phase velocity curves (fundamental plus upto 3th higher mode) on the right are shown.
The researcher in [33] has followed the (p,τ) stack method proposed by a researcher to obtain the surface wave phase velocity dispersion from an array of seismic traces. Here p and τ stand for the slowness and intercept time, respectively [34]. The Fourier spectrum of a seismic signal observed at distance r n is given as follows.
( , ) ( ) eq. (1) Here is the amplitude, is the phase, is the angular frequency and = √−1 is the complex number. The following relation defines the ( , ) stack of N traces at different distances from the same source.
The latter is used to correct the observed surface wave amplitudes for the geometrical spreading and to remove the source phase. The maxima of |F(p,ω)| is searched to obtain the possible dispersion curves.

Results and Discussion:
We measured a total of 33 resistivity profiles. The individual profiles close to linear are mostly 155-m long, but some of them are shorher due to setting the electrode spacing to a value smaller than 5-m. The total profile length is then about 4-km (see Fig. 6). For the same profile setting we measured surface waves for the MASW analysis. We will not present all our measurements herein, but will share a map summarizing our findings for the underground void structure beneath the Nevşehir Castle region.
We use the rainbow colors to present our tomography results where red, green and blue colors stand for high, medium and low values, respectively. Where appropriate black filled circles are used to indicate the locations of isolated underground structures interpreted as voids at certain depth. The physical size of these voids (i.e. width, length and depth) will be interpreted using the electrical resistivity imaging sections. Because of the low resolution of the seismic measurements due to the lack of high frequency energry in the observed fundamental mode Rayleigh surface waves, information regarding particulary depth to the void may be misleading on the seismic sections. The seismic depth information is better confined by the use of higher modes, but we were not able to employ these modes in our analysis due to highly scattering nature of underground structures beneath the studied area.
In the presence of an underground void, the corresponding electrical resistivity values are high (i.e. red color) while the seismic velocity values are low (i.e. blue color). In this way, the results of electrical resistivity and seismic velocity should support each other. On the other hand, if the electrical resistivity and seismic velocity values are both high (i.e. red color), then a dry, strong, not deformed, most likely magmatic origin rock mass beneath the survey line at some depth is expected. On contrary, if the electrical resistivity and seismic velocity values are both low (i.e. blue color), then, in this case, a wet, weak with cracks, most likely sedimentary origin rock mass is expected to be the cause of the observed anomaly.
In Fig. 9, we show an example of fundamental-mode Rayleigh surface wave phase velocities extracted by using the "Computer Programs in Seismology" by [33]. The left panel shows the phase velocity stack values (see eq. 2) obtained from an array of seismic traces (currently N=24). The color coding corresponds to the stack values where the red color gives the highest stack value and the black symbols indicate the peaks in the stack described over the phase velocityperiod grid. The right panel gives the result of 1-D inversion of the observed Rayleigh fundamental-mode phase velocity curve chosen from the highest value stack in the left panel. The inverted 1-D shear-wave velocity-depth profile (small inset on the right) has a depth range down to 42 m where the shear velocity is ~1.0 km/s and indicates a velocity reversal around 27-m depth. The phase velocity dispersion curves (fundamental plus 9 higher modes) predicted from the inverted model overlay the phase velocity stack values in the left panel from which it is noted that the higher modes are difficult to identify. The same procedure is applied to all other sets of seismic traces in the current work.

Figure 9:
On the left is shown phase velocity stack values from eq. (2) and on the right is the theoretical fit to the observed fundamental-mode phase velocities along with the inverted shear-wave velocity-depth profile.

Profile Eski Gore Yolu (sub-profiles 11 and 12)
The profile length is 310-m. The resistivity measurements each with 32 electrodes are taken in two parts. The upper panel in Fig. 10 shows the electrical resistivity cross sections. The three high resistivity depths (two in the first cross section on the left and one in the second cross section on the right) are interpreted as underground voids at depths deeper than 10-m and shallower than 25-m. The first and the second voids are around 38-m and 112-m distances in the first profile and the third void is around 62-m in the second profile.
The seismic inversion results are shown in the lower part of Fig. 10. The seismic profile length is 368-m. Eight 24-geophone spreads cover this profile and are interpreted together for a 2-D seismic cross section. The low velocity zones depicted by the blue color depths in the seismic cross section well correlate with the high resistivity zones depicted by the red color depths in the resistivity cross sections. The three arrows and the black filled circles are used to show this correlation. The bedrock depths in the seismic cross section are depicted by the red color zones. It appreas that the early settlers avoided digging the voids in the bedrock while effectively utilizing the relatively low velocity (soft) rock at shallower depths. This situation is found valid for all other profiles measured. Note that these void structures in the seismic cross section appear deeper than they appear in the resistivity cross sections. This is due to the fact that in the fundamental-mode dispersion curve the deeper structures are modeled by only longer periods with relatively lower resolution. The higher mode dispersion curves would better model the deeper structures, but we were not able to include them in the inversion. We think that the anomaly depths given by the resistivity cross sections are currently robust. In Fig. 10

Profile Hapishane Yolu (sub-profiles 13 and 14)
This profile is almost parallel to the previous profile (Eski Gore Yolu; see Fig. 6), but there is ~30-m difference in altitude. The total length of this profile is 310-m and the resistivity measurements are taken again in two parts. In Fig. 11, the correponding two resistivity cross sections in the upper panel and the seismic velocity cross section in the lower panel are shown. The two void structures on the left and another one on the right observed in the resistivity cross sections for the profile "Eski Gore Yolu" above are similarly observed for these profiles as well. These void structures take place in the depth range 8-20 m. The first and the second voids are positioned around 30-m and 105-m in the first profile and the third void is around 38-m in the second profile. The void structures appear to gain altitude paralel to the topography going upward from the profile "Hapishane Yolu" to the profile "Eski Gore Yolu".
The length of the seismic profile is again 368-m. Eight shot gathers are combined to interpret the underground seismic structure in 2-D. The three arrows and the black filled circles are used to indicate the correlation between the resistivity and seismic results. The high resistivity values (red color) on the resistivity cross sections correspond to the depths of low-velocity zones (blue color) on the seismic cross section. The bedrock depicted by the varying tones of the red color is relatively shallower. The relation between the void structures described in Figs. 10 and 11 is better explained in Fig. 12 where the resistivity cross sections corresponding to the profiles "Eski Gore Yolu" and "Hapishane Yolu" are shown together. Even though these profiles are well separated from each other in the vertical and horizontal directions, it is easy to observe the extension of the void structure beneath the profile "Eski Gore Yolu" to continue under the profile "Hapishane Yolu". If our interpretation is correct regarding the findings in Fig. 15, then properly positioned other profiles parallel to the profiles in Figs. 13 and 14 will provide extra information for the extension of the tunnels going up at the site of the Nevşehir Castle's twin and down to the plain terrain (see Fig. 6).

Figure 12:
The resistivity cross sections for the profiles "Eski Gore Yolu (sub-profiles 11 and 12)" and "Hapishane Yolu (sub-profiles 13 and 14)" are shown together. The cross sections are viewed from the flat terrain looking toward the Nevşehir Castle's twin.

Profile Bas Cesme (sub-profiles 8, 9 and 10)
This profile is the longest one with 465-m length. The resistivity measurements were taken in three parts. The three resistivity cross sections shown in Fig.  13 (upper panel) indicate three void structures separately evident in each cross section. The void depths are shallower than 20-m and deeper than 8-m. The first void in the first profile is at ~85-m distance, the second one is at ~37-m in the second profile and the third one is at ~100-m in the third profile.
The profile length for the corresponding seismic cross section is 460-m. The 2-D underground seismic structure is interpreted using ten shot gathers. The high resistivity zones depicted by the red color depths in the resistivity cross sections well correlate with the low velocity zones depicted by the blue color depths in the seismic cross section. This correlation is indicated by the three arrows and the black filled circles. The bedrock depicted by the varying tones of the red color lays almost horizontally. Other parallel profiles on the up-hill and down-hill side of the profile "Bas Cesme" (particularly the sub-profiles 8 and 9) may help us further track the void structure beneath the surface (see Fig. 6). However, due to the restrictions created by ruins after the demolition works we were not able to measure extra profiles in this part of the Nevşehir Castle region.

Profile Kaya Cami Ust Sokak (sub-profiles 4 and 5)
The total length of this profile is 310-m. The resistivity measurements were taken in two parts for which the corresponding resistivity cross sections are shown in the upper panel of Fig. 14

Several Profiles Interpreted Together
Close inspection of the anomalies on the resistivity and seismic cross sections show that there exist systematic void structure carved in the underground within certain depth range, which is consistent with the void structure revealed by the actual clean-up work shown in Fig. 1 (lower panel). Those thick lines in Fig. 15 almost parallel to the topography of the cone shaped Nevşehir Castle region indicate the interrelationship of the void system from one cross section to another. Fig. 15 shows the resistivity cross sections of several profiles together for a better visual inspection. These profiles take place on the northeast flank of the Nevşehir Castle region, i.e. the sub-profile 10 in the profile "Bas Cesme" and the profiles "Kaya Cami Ust Sokak (i.e. sub-profiles 4 and 5)", "Kaya Cami (i.e. sub-profiles 1, 2 and 3) " and "Vaiz Osman (i.e. sub-profiles 6 and 7)" (see Fig. 6).
Due to the limitations in terms of time and resources we had to set the resistivity electrode spacing to 5-m, which is relatively larger than needed to resolve the fine details of the void structure of such a large area. Those voids 1-3 m wide are weakly resolved, but those ones 5-6 m wide or more are resolved well in the observed resistivity and seismic data. The clean-up work shown in Fig. 1 (lower panel) indicates that the underground settlements beneath the Nevşehir Castle region has a multiple-story structure going upward from the flat terrain to the castle on top where each level is represented by many entrances distributed on a perimeter with more or less same altitude. The three cross sections (i.e. sub-profiles 10, 1 and 4) are shown in the upper left of Fig. 15 where the three high resistivity depths (red color) are void structures beneath each profile. The thick black line (approximately 125-m long) shows the interrelationship between these voids at three different locations in a 3-D sense.
A similar interpretation is given for the high resistivity depths under the sub-profiles 6, 2 and 4 (upper right). Another thick black line approximately 125m long is used to show the interrelationship between the void structure on this side of the studied region. Note that the resistivity anomaly under the subprofile 6 appears relatively weakly developed. The void under this location is partially filled with sand, which is an information confirmed with Nevşehir City officials. It was said that after collapse in the road the resulting trough was filled with a few truck loaded sand (i.e. a few tons). The 3-D interpretation is extended for the void structure under the sub-profiles 7, 3 and 5 (lower panel) where again the thick black line (approximately 129-m long) represents the respective interrelationship between those voids on different levels of the underground settlements. The void structures under the sub-profiles 7 and 5 are well developed compared to the one under the sub-profile 3 that appears relatively smaller. Perhaps the void gets narrower in this location or it is partially filled with other earth material along with or without a downfall. Extra profiling especially within galleries with narrower spacing of electrodes and geophones are needed to locate the inner (or completely hidden) voids more correctly in three-dimensions. Other geophysical methods such gravity, magnetics and ground penetrating radar may help for a more complete geophysical investigation of the region. The important question is that how these voids continue after reaching the flat terrain on which the city of Nevşehir occupies most of the area.

Conclusion:
All the measurements taken after the "Urban Transformation" and before the "Clean-up Work" (see Fig. 1) are interpreted and the results are summarized in Fig. 16 where the profile numbers in red color refer to those profile numbers given in Fig. 6. Interrelationships between these voids with a depth range approximately 8-20 m are indicated using the green lines. Note that part of the overlying material (man made and soil) once visually hiding the voids from the view is removed so that entries to those voids are currently visible. The various measurement profiles encircling the study region cross over the underground void structures radially emanating from the Nevşehir Castle's top going to the plain terrain. The underground city under consideration has a multiple-story construction with many entries and has much more detailed interior structure with connecting galleries than reflected herein with a limited number of observations. The following list provides some recommendations for the future work.
• Underground cities such as the one currently studied herein belong to the world heritage. Therefore, every effort should be made to unearth these precious structures to carefully convey them to the next generations.
• Geophysical methods should be applied especially in an integrated manner so that the underground structures are discovered with great precision while protecting them from undesirable harms during the discovery.

•
Collaboration with other disciplines such as geology and archaeology should increase the effectiveness of the future work that should be planned for the Nevşehir Castle region and the surrounding.