MODELLING SPELEOGENESIS IN SOLUBLE ROCKS: A CASE STUDY FROM THE PERMIAN ZECHSTEIN SEQUENCES EXPOSED ALONG THE SOUTHERN HARZ MOUNTAINS AND THE KYFFHÄUSER HILLS, GERMANY MODELIRANJE SPELEOGENEZE V TOPNIH KAMNINAH: PRIMER ZECHSTEINSKIH KAMNIN V JUŽNEM DELU GOROVJA HARZ IN V HRIBOVJU KYFFHÄUSER Georg KAUFMANN 1, * & Douchko ROMANOV 2 Abstract UDC 551.435.8(430) Georg Kaufmann & Douchko Romanov: Modelling speleo- genesis in soluble rocks: A case study from the Permian Zech- stein sequences exposed along the southern Harz Mountains and the Kyffhäuser Hills, German Soluble rocks such as limestone, dolomite, gypsum, anhydrite, and salt can be dissolved by water flowing through voids in the rocks. The removal of the dissolved material from fissures and bedding partings by physical and/or chemical dissolution enlarges the permeability of the soluble rocks within geologi- cally short periods of time, ranging from 100,000 years down to decades. This geologically short evolution time of voids in soluble rocks poses a substantial risk of mechanical instability of the enlarged voids, and possible surface deformation, when enlarged voids start to collapse. We describe karst and cave features in the rock sequence ex- posed along the southern part of the Harz Mountains and the Kyffhäuser Hills in Germany, where limestone/dolomite and anhydrite/gypsum are exposed along a kilometer-wide strip following the foothills of the Harz Mountains. The rocks have been deposited during the Permian Zechstein period, buried, and exposed later through tectonic uplift. The exposed part of this soluble sequence is dominated by karst features. But there are also substantial cave voids deeper in the rock, with no ob- vious entrance to the surface, which have been discovered by chance through mining activities. Often, the sub-surface void evolution is closely linked to surface deformation, creating col- lapse sinkholes and subsidence. In the city of Bad Frankenhau- sen at the foothills of the Kyffhäuser Hills, the evolution of sub- Izvleček UDK 551.435.8(430) Georg Kaufmann & Douchko Romanov: Modeliranje spe- leogeneze v topnih kamninah: primer zechsteinskih kamnin v južnem delu gorovja Harz in v hribovju Kyffhäuser Voda, ki teče skozi pore in razpoke v topnih kamninah, kot so apnenec, dolomit, sadra, anhidrit in sol, raztaplja stene prevod- nih poti. Učinkovito odnašanje raztopljene snovi lahko močno poveča hidravlično prevodnost vodonosnika v geološko krat- kem času, ki je v izjemnih primerih dolg vsega nekaj desetletij. Hiter razvoj prevotljenosti lahko povzroči mehansko nestabil- nost nastalih votlin in posledično ugrezanje površja. Tak prim- er najdemo v kilometer širokem pasu apnenca, dolomita, sadre in anhidrita, ki se razteza vzdolž južnega dela pogorja Harz in hribovja Kyffhäuser. Sedimenti, ki so se odložili v permskem Zechsteinskem morju, so bili kasneje globoko pokopani in ponovno tektonsko izdani na površje. Izdanki teh kamnin so izrazito kraško preoblikovani, prevotljenost pa je velika tudi v globlje pokopanih kamninah, ki nimajo očitne povezave s površjem in so bile odkrite pri rudarskih delih. Velikokrat se razvoj votlin pod površjem izrazi tudi na površju, kjer nastaja- jo udornice in grezi. V mestu Bad Frankenhausen ob vznožju hribovja Kyffhauser se je zaradi razvoja votlin pod površjem nagnil cerkveni zvonik. V članku raziskujemo razvoj kraških si- stemov v apnencu in anhidritu z numeričnimi modeli. V mod- elu upoštevamo tok, raztapljanje in prenos snovi v kamninskem masivu, sestavljenem iz topnih in netopnih kamnin, kjer razta- pljanje poteka v apnencu in anhidritu. Ključne besede: topne kamnine, kras, razvoj jam v rudnikih, udorne doline, numerično modeliranje. 1,* Georg Kaufmann, Institute of Geological Sciences, Geophysics Section, Freie Universität Berlin, Malteserstr. 74-100, Haus D, 12249 Berlin, Germany, e-mail: georg.kaufmann@fu-berlin.de 2 Douchko Romanov, Institute of Geological Sciences, Geophysics Section, Freie Universität Berlin, Malteserstr. 74-100, Haus D, 12249 Berlin, Germany, e-mail: douchko.romanov@fu-berlin.de Received/Prejeto: 25.02.2019 DOI: https://doi.org/10.3986/ac.v48i2.7282 ACTA CARSOLOGICA 48/2, 173-197, POSTOJNA 2019 COBISS: 1.01 GEORG KAUFMANN & DOUCHKO ROMANOV INTRODUCTION Soluble rocks from the Permian Zechstein period are present underneath large parts of Northern Germany. These Permian Zechstein deposits comprise sequences of soluble rocks, with limestone along the base of a se- quence, followed by anhydrite and salt, and often topped with insoluble marls and/or schists. While in most parts of Northern Germany these Permian Zechstein sequenc- es can be found in several kilometer depth, in some parts the soluble rocks have been moved closer to the surface by salt diapirism and tectonic uplift. Along the southern Harz Mountains in Germany, the Permian Zechstein sequences are exposed as a result of the tectonic uplift of the Harz along a narrow kilome- ter-wide strip, with widespread karstification within the outcropping soluble rocks. The uplifted soluble rocks dip gently towards the west / south-west, and karst features such as collapse sinkholes and resurgences are also mark- ing the buried part of the dipping rocks. In the past, the base of the Zechstein sequence, a copper-shale ore (Kup- ferschiefer), has been mined extensively, with the discov- ery of numerous karst features such as a special type of large mine caves (e.g., Brust & Nash 2013; De Waele et al. 2013), termed Schlotten in german, in the normally impermeable anhydrite layers. Besides these often spectacular karst features, the dissolution of the soluble rocks (limestone/dolomite, an- hydrite/gypsum, salt) is responsible for subsidence of the surface and resulting collapse sinkholes, once the over- burden covering the enlarged voids evolved in the soluble rocks becomes too weak. A spectacular example of the ongoing dissolution in the Permian Zechstein sequence underneath the city of Bad Frankenhausen (Kyffhäuser Hills) is the inclined church tower of the Oberkirche, which has been build on a large fault with Permian Zech- stein rocks in shallow depth. The titling of this church tower, with an inclination amounting to 4.9 o , is caused by dissolution of the soluble rocks in the sub-surface, a still active process in the region. We are interested in the solution process in the solu- ble rocks, especially the limestone and anhydrite layers in the shallow sub-surface of the Permian Zechstein rocks. As the timescales for dissolving limestone and anhydrite differ by more than one order of magnitude, interesting interactions between these types of soluble rock occur. We will study, by numerical means, a simplified setup of the inclined Permian Zechstein rocks, with dissolution in the limestone, dissolution in the anhydrite with associ- ated precipitation of gypsum, and a regional water flow providing aggressive solution from below (hypogene set- tings). We organize this paper as follows: In the second section, we introduce the broader setup of the region under consideration, discuss the sur- face outcrop of the Permian Zechstein sequences with its soluble rock types carbonates, sulfates, and halites, and their continuation in the sub-surface further south. We will point out the hypogenic origin of many karst features in the region. In the third section, we discuss three caves typical for the region, all developed in the anhydrite layers of the Permian Zechstein sequences, and their relation to sur- face features such as sinks, resurgences, and sinkholes, the latter related to significant surface deformation and thus a threat to infrastructure. In the fourth section, we present the chemical reac- tions for the interaction of water and its dissolved agents with the soluble rock present in the working area. W e dis- cuss chemical equilibria, the way towards these chemical equilibria, and present time scales for dissolving rocks types such as the carbonates and the sulfates of the Perm- ian Zechstein sequences. W e finally discuss the combined dissolution of anhydrite and the precipitation of gypsum, usually responsible for the sealing of the anhydrite as aquitard. In the fifth section, we use our chemical derivation and embed it into a three-dimensional flow, transport and reaction model, describing a simplified version of the Permian Zechstein rocks (only one sequence). We will see what is needed to generate Schlotten-type mine caves in the anhydrite rock. Our key questions in this paper are: 1. Can we explain the evolution of the Schlotten-type surface voids is responsible for the tilting of the church tower of the Oberkirche. We explore the evolution of such a karst system composed of limestone and anhydrite by numerical means, describing flow and transport in a rock mass composed of soluble and insoluble rock sequences, with limestone and anhydrite responsible for the evolution of secondary porosity. Key words: soluble rocks, karst, mine-cave evolution, collapse sinkholes, numerical modelling. ACTA CARSOLOGICA 48/2 – 2019 174 MODELLING SPELEOGENESIS IN SOLUBLE ROCKS: A CASE STUDY FROM THE PERMIAN ZECHSTEIN SEQUENCES EXPOSED ALONG THE SOUTHERN HARZ MOUNTAINS AND THE KYFFHÄUSER HILLS, GERMANY mine caves with a large-scale karst evolution model, accounting for dissolution in a combined carbonate- sulfate setting? 2. Can we explain the surface deformation observed in Bad Frankenhausen with numerical models of large- scale dissolution in the underlying sulfate layers? 3. Can we estimate the role of aggressive mineral-rich brine water arriving from below through fault systems (hypogene setting)? GEOLOGY In this section, we introduce the geological setting of the working area and discuss the lithology in the broader context of the soluble rocks exposed along the Southern Harz Mountains. HARZ MOUNTAINS The Harz Mountains in Northern Germany (Fig. 1) are mainly composed of old magmatic and metamorphic rocks (e.g., Hohl 1985; Schönenberg & Neugebauer 1987), and are surrounded by soluble rocks, mainly Cretaceous limestone in the North and Permian anhydrite and lime- stone/dolomite in the South. The mountain ridge, about 110 km long and 30-40 km wide, with its highest peak the Brocken (1141 m), has been uplifted during the V ariscian mountain building period (700-500 Ma B.P .), and was al- most denuded at the beginning of the Permian Period, when the Permian Zechstein ocean transgressed several times over this area (300-250 Ma B.P.), which was then located on the equator. The repeated deposition from the Zechstein ocean resulted in up to seven depositional cy- cles (z1/zW-Werra, z2/zS-Staßfurt, z3/zL-Leine, z4/zA- Aller, z5/zO-Ohre, z6/zFr-Friedland, z7/zFu-Fulda), all starting with carbonates (limestone/dolomite), sulfates (anhydrite, gypsum), halides (halite, sylvite), and finally shales. The bottom of the Permian Zechstein series is the Kupferschiefer, which has been mined in the past exten- sively. After the Zechstein ocean finally disappeared, ero- sion deposited the late Permian and Triassic Buntsand- stein rocks, topped by Muschelkalk. During the uplift of the European Alps (100 Ma B.P .), the Harz region got up- lifted again, but erosion lowered the surface since then. A final uplift phase in the Tertiary caused an asymptotic uplift of the Harz Mountains, with the northeastern part forming a steep cliff and the southwestern part experi- encing less uplift. In the wake of this tilting, the Perm- ian Zechstein rocks, already buried at greater depth, were uplifted to the surface along a narrow stripe following the Fig. 1: Topographical map of the Harz Mountains (with the Brocken the highest peak) and the Kyffhäu- ser Mountains (between Kelbra and Bad Frankenhausen). Cities are shown as grey dots, landmarks as grey triangles. The grey-hashed area marks the outcrop of the solu- ble Zechstein rocks. In the inset, the location of the map within Germa- ny is shown as red rectangle. ACTA CARSOLOGICA 48/2 – 2019 175 GEORG KAUFMANN & DOUCHKO ROMANOV Southern Harz. This area nowadays is intensively karsti- fied. KYFFHÄUSER HILLS The tectonic evolution of the Kyffhäuser Hills (Fig. 2), a 6x13 km outcrop south of the Harz Mountains, is inti- mately linked to the Harz Mountains (e.g., Wunderlich 2014). The Kyffhäuser Hills, with its highest peak the Kulpenberg (473 m), experienced a similar asymmetric tilting during Tertiary times, creating a steep northern edge and exposing Permian Zechstein rocks along the southern rim. Additionally, two major faults bound the Kyffhäuser Hills, the northern (NKF) and the southern (SKF) Kyffhäuser Fault. Along the SKF, the Permian Zechstein rocks have been displaced vertically by several hundreds of meters. GEOLOGICAL CROSS SECTIONS Next, we discuss the lithological structure of the en- tire region. In Fig. 2, three profiles are marked with red lines. Geological cross sections of these three profiles are shown in Fig. 3. The first section ( A-A ’-A ’’, Schriel & Bülow 1926b) traverses the southern Harz Moutains close to the city of Sangerhausen. Here, the Kupferschiefer has been mined extensively, and the lithologies are well known from min- ing. In Fig. 3 (top), the three letters mark mining shafts, with A the location of the Röhrigschacht, A’ the location of the Thomas-Mün(t)zer-Schacht, and A’’ the location of the Bernhard-Koenen-Schacht. Starting from the north, the Permian Zechstein sequences z1/zW (Werra) and z2/ zS (Stassfurt) are exposed along the surface, caused by the tilting of the Harz. The soluble sequences then gently dip towards the south, and are covered by younger de- posits, amongst other the Buntsandstein sequence. While the salts are absent in the sequences for the depth range from surface to about -200 m below sea level, they can be found further south (of Thomas-Münzer-Schacht), where the Permian Zechstein rocks are buried deeper and the salt layers have not been removed by circulating water from the surface, which is devoid of dissolved car- bonates, sulfates, and salts. This cross section is represen- tative for a location along the southern Harz Mountains, in which the Permian Zechstein layers are more or less undisturbed by further faulting. In the second section in Fig. 3 (middle, B-B’, Wun- derlich 2014), the three first Permian Zechstein sequenc- es are present (z1/zW-z3/zL). While in the western part close to the Kyffhäuser Hills the salts have been dissolved by circulating water, in the eastern part, beyond the Northern Kyffhäuser Fault (NKF), the salts are still pres- ent. The lithologies are confirmed by the deep boreholes Ud17/56 and Ud119/64. The third section in Fig. 3 (bottom, C-C’, Schriel & Bülow 1926a,b; Wadas et al. 2016) crosses the south- ern Kyffhäuser Hills towards the south. Here, north of the Southern Kyffhäuser Fault (SKF) the Permian Zech- stein sequences z1/zW and z2/zS are exposed along the surface. South of the SKF, these sequence can be found around 300 m deeper, offset by the Southern Kyffhäuser fault, and here the salt in the z2/zS-Leine Series is still present. Thus, an interesting situation arises: the exposed soluble rocks north of the SKF are intensely karstified by meteoric waters, with numerous dolines, sinkholes, and caves, and along the SKF brackish water from deeper lay- ers ascends and increases the calcium equilibrium con- centration in the layers above due to the presence of dis- solved salt ions. FROM GEOLOGY TO KARST … We have seen that the entire region (Southern Harz Mountains, Kyffhäuser Hills) is characterized by solu- ble rocks of the Permian Zechstein Period, which out- crop along a narrow strip due to the tectonic uplift of the Southern Harz Mountains and the Kyffhäuser Hills, and then continue into the sub-surface with a shallow dip an- gle towards the south-west. Faulting has altered this sim- ple sequence at some locations, with significant impact to the karstification, as we will see in the next section. Fig. 2: Topographical map of the Kyffhäuser Mountains and its vicinity. The blue area is the Kelbra Reservoir, the black dashed lines the northern and southern Kyffhäuser fault (NKF, SKF) and the Kelbra Fault (KF), and the red lines with markers the profiles shown in Fig. 3. The three caves marked with bold letters and red polygons are: ES: Elisabethschächter Schlotte, BC: Barbarossa- Cave, NC: Numburg-Cave. ACTA CARSOLOGICA 48/2 – 2019 176 MODELLING SPELEOGENESIS IN SOLUBLE ROCKS: A CASE STUDY FROM THE PERMIAN ZECHSTEIN SEQUENCES EXPOSED ALONG THE SOUTHERN HARZ MOUNTAINS AND THE KYFFHÄUSER HILLS, GERMANY KARST In this section, we present evidence for the widespread karstification of the region, with its special mine caves, the Schlotten, entrance-less voids in the anhydrite lay- ers of the Permian Zechstein rocks, often exposed by mining operations. We then link the knowledge about sub-surface karst features to the situation in the city Bad Frankenhausen, where the tilted church tower of the Oberkirche is a spectacular sign of solution processes in the sub-surface. CAVES In the soluble Zechstein rocks, numerous caves can be found. A special type of caves are the Schlotten, mine caves in the anhydrite rock just above the limestone or dolomite layer. These Schlotten classify as hypogenic caves (e.g., Klimchouk 2011; Kempe 2008; Bauer & Elste 2016; Kempe et al. 2017), with water supply responsible for dissolution arriving through the underlying lime- stone/dolomite layers. Fig. 3: Geological cross sections for the three profiles shown in Fig. 2. Faults are shown as black lines and marked (northern and southern Kyffhäuser fault, NKF, SKF), boreholes as white bars with black outlines. The letters indicate positions marked in Fig. 2. ACTA CARSOLOGICA 48/2 – 2019 177 GEORG KAUFMANN & DOUCHKO ROMANOV These caves are often found by mining operations and are mainly located below the local base level, but ac- cessible due to the pumping in the mine galleries. The german term Schlotte is a mining term describing natural cavities found during the mining operations, often in the anhydrite layers of the Permian Zechstein sequences. The Schlotten have often been used by the miners as a natural drainage route for mine water, and to store waste rock. Next, we discuss three prominent caves of this type, which are characteristic for the karstification of the re- gion. Elisabethschächter Schlotte The Elisabethschächter Schlotte is a large mine cave in the Werraanhydrit (z1/zW) (Fig. 4), found through mining operations during mining of the Kupferschiefer already in 1755, it has no natural entrance (e.g., Völker & Völker 1982; Brust 2005; Kupetz & Brust 2008; Kupetz 2008; Spilker 2016; Brust & Graf 2016). It is located just south of the exposed soluble rocks along the southern Harz rim (see marker in Fig. 2), in a classical undisturbed situation of the soluble Permian Zechstein rocks. The cave is around 400 m long, with a succession of large rooms up to 14 m high. The cave follows the in- clined bed of the Zechsteinkalk (z1/zW) layer, thus the rooms are dipping towards the south. Cave rooms, how- ever, are developed in the Werraanhydrit (z1/zW). The larger cave rooms, often cupola-shaped, are connected by small crawlways. Most of the floor of the cavity is a fine-grained sediment, made up by the insoluble parts of the anhydrite. In the northern part, the cave ends in a large breakdown room, the Chaosdom. Here, a connec- tion to a sinkhole on the surface is evident, as numerous surface materials such as leaves, branches, and even arti- ficial building material slowly migrate into the cave. The distance from the top parts of the Chaosdom to the sur- face, however, is still 44 m. The southern part of the cave is its lowermost section and the last room Kuppeldom is completely filled with fine-grained sediments. While the cave extends vertically between 174-221 m N.N., the local base level Grenzbach is higher (266 m N.N.). The cave interior is characterised by debris of gypsum mate- rial. This material is a result from the anhydrite-to-gyp- sum conversion in the cave rooms: The high humidity dissolves anhydrite from cave roof and walls, but gypsum is (almost) immediately precipitated. The increase in volume during the gypsum precipitation results in loose layers, which then easily detach under the gravitational acceleration imposed on them. With the schematical cross section (Fig. 4), we de- scribe the current hypothesis for the evolution of the cave: Aggressive water from the surface sinks along the permeable fractured limestone aquifer (in this case the Zechsteinkalk) and reaches calcium saturation with re- spect to the limestone. As the water table (before mining started) has been way above the cave location, the cavities evolved under deep phreatic conditions. Here, along the upper interface of the Zechsteinkalk, the water is in con- tact with Werraanhydrit, thus the water can still dissolve anhydrite along the interface. The widespread dissolution of the Werraanhydrit can be traced in the mine galleries along the limestone-anhydrite interface, as here the in- soluble residuals from the Werraanhydrit accumulate as fine-grained sediments. However, the evolution of larger cupola-shaped rooms cannot be explained so easily, as they are not so widespread. Here, additional conditions need to be met, such as structural controls (e.g., faults) or preferential flow in the Zechstein layer causing an in- homogeneous distribution of water along the interface. Note that the cave follows the steep interface Zechsteink- alk-Werraanhydrit. Barbarossa-Cave The Barbarossahöhle is another mine cave of the typus Schlotte, also located in the Werraanhydrit (z1/zW) just above the Zechsteinkalk (z1/zW) layer (Fig. 4). The Bar- barossa-Cave is located south of the Kyffhäuser Hills in the Permian Zechstein rocks, with the southern Kyffhäu- ser fault (SKF, for location see Fig. 2) nearby. This large cave has been found in 1865 in search for new locations to mine the Kupferschiefer at the base of the Permian Zechstein rocks (e.g., Brust 2005; Kupetz 2005; Adler & Mertmann  2012; Kupetz 2014). Soon af- ter the discovery, the mining operations were abandoned and the cave opened as a show cave. The cave is around 800 m long, mainly horizontal and is characterised by a succession of larger rooms, often with wide cupola-shaped roofs. The groundwater table is exposed at various places, and in contact with the local base level. The spring Pfannspring about 700 m north of the cave provides surface water, which sinks after a short distance and reappears in the Barbarossa-Cave. Flow ve- locities, however, are very small (Adler 2013). High humidity causes conversion of the anhydrite to gypsum along the roof areas of the chambers, the vol- ume increase causes the layers of converted gypsum to shear off and fall down, once they become too heavy. This process is responsible for room enlargement. Parts of the large chambers are already in a state of breakdown, e.g., the Olymp, which reaches close to the surface, were the collapse sinkhole Teufelsgrube is located. The Barbarossa-Cave does not follow the inclined bed of the Zechsteinkalk (as in the case of the Elisabeth- schächter Schlotte), instead it develops horizontally to- wards the local base level, the Wipper valley. The cave has formed under shallow phreatic conditions, but still its ACTA CARSOLOGICA 48/2 – 2019 178 MODELLING SPELEOGENESIS IN SOLUBLE ROCKS: A CASE STUDY FROM THE PERMIAN ZECHSTEIN SEQUENCES EXPOSED ALONG THE SOUTHERN HARZ MOUNTAINS AND THE KYFFHÄUSER HILLS, GERMANY Fig. 4: Cave maps (left) and schematical cross sections (right) of three caves. The cave maps are shown in red, superimposed on the topogra- phy (color coded). For legend of geological cross sections see Fig. 3. Cave cross sections are shown in white, hypothetic water table as dashed black line, and faults as solid black lines. Top: Elisabethschächter Schlotte, Middle: Barbarossa-Cave (dashed circle, location of collapse sinkhole), Bottom: Numburg-Cave (solid black outline-mapped parts, dotted black outline-hypothetical part). For legend see Fig. 3. ACTA CARSOLOGICA 48/2 – 2019 179 GEORG KAUFMANN & DOUCHKO ROMANOV origin is hypogenic water from the Zechsteinkalk layer below the Werraanhydrit. Numburg-Cave The Numburghöhle is located east of the Kyffhäuser Hills along the southern shore of the Kelbra Reservoir (see Fig. 2 for location). The cave, which is developed in the Werraanhydrit (z1/zW), has been mapped to 1750 m length, and is mainly developed below the (current) wa- ter table (Völker and Völker 1991, 2014). The entrance of the cave has long been known as a cave spring, with brackish water emerging from a cave room, ending in a syphon. Starting in 1987, the exten- sion of the Kupferschiefer-mining (Thomas-Mün(t) zer-Schacht, Bernhard-Koenen-Schacht) towards the west with its large-scale pumping was hampered by a substantial increase of leakage into the mining galleries. The groundwater withdrawal tapped the region around the Kelbra Reservoir, and the source of the Numburg- Cave became a sink, sucking around 10 m 3 of water per minute from the reservoir. The water level in the cave started to drop, and in August 1988 the phreatic parts of the cave became accessible, revealing a large cave (Fig. 4). As the inflow of water into the mining galleries become too large to be managed, the western part of the mine was sealed off by closing galleries with sealing dams and the water level in the Numburg-Cave and its surround- ing started to rise again. In April 1989, the cave become flooded again up to its initial height. The Numburg-Cave is a succession of large rooms, and it ends abruptly along the Kelbra fault. Here, in the last cave room, sandstone from the Carbon Period is ex- posed along a fault zone, and from this rock, water in- filtrates the cave from several fissures in the sandstone. In the cave, brackish water emerges through two shafts, which have been explored to around 22 m depth. This brackish water ascends along the Kelbra fault zone from the Zechstein deposits in the north, which are located in greater depth, with the salt still present. Thus, the Numburg-Cave has at least two sources of water able to dissolve the Werraanhydrit, the aggressive water from the carboniferous sandstone, and the brack- ish water rising along the fault, the latter increasing the solubility of gypsum/anhydrite because of the presence of the dissolved salt. Bad Frankenhausen We now move on from the known caves and cavities in the soluble Permian rocks of the Southern Harz and Kyffhäuser region to the city of Bad Frankenhausen, which is located right on the Southern Kyffhäuser Fault (SKF). In Bad Frankenhausen, the church Oberkirche has been build in 1382, but around 1640 the tower started tilting towards the east. Nowadays, the leaning church tower has an inclination of 4.9 o (Fig. 5). The reason for this inclination can be found in the sub-surface: As al- ready mentioned, the Southern Kyffhäuser Fault runs just south of the church, and here the soluble Permian Zechstein rocks are offset vertically by about 200 m. Sev- eral sinkholes, clustered around the Southern Kyffhäuser Fault, can be found, with the largest ones the Quellgrund, an old sinkhole in Bad Frankenhausen known since me- dieval times (998, and currently the spa gardens), and the Äbtissinnengrube, an elongated sinkhole 160 by 120 m in size, known since the 16 th century. On the western site of the latter sinkhole, in 2009 a new sinkhole appeared, 20 m in diameter and 13 m deep, indicating the active dissolution processes (e.g., Brust 2010). Since 1997/1999, the area around the church tower has been surveyed by precise levelling (Scheffler & Mar- tienßen 2013). Two points from the 8 node levelling net- work show significant horizontal movements of around 74-90  mm and vertical subsidences of around 100- 160 mm within the period 1997/1999-2012 (Fig. 6). The Fig. 5: Leaning church tower of the Oberkirche in Bad Franken- hausen (Photo: G. Kaufmann). ACTA CARSOLOGICA 48/2 – 2019 180 MODELLING SPELEOGENESIS IN SOLUBLE ROCKS: A CASE STUDY FROM THE PERMIAN ZECHSTEIN SEQUENCES EXPOSED ALONG THE SOUTHERN HARZ MOUNTAINS AND THE KYFFHÄUSER HILLS, GERMANY deformation is, however, not uniform in time. Since 2014, Kersten et al. (2017) conduct both geodetic (120 points) and gravimetric (13 points) measurements in Bad Fran- kenhausen. Annual local subsidence values measured during this campaign amount to 4-10 mm, however, is partly related to recent construction work, thus separa- tion of natural and anthropogenic effects is still difficult. In 2013 and 2014, several boreholes have been lowered in vicinity of the church tower by the TLUG (Thüringer Landesanstalt für Umwelt und Geologie) and in parts by the LIAG (Leibniz-Institute of Applied Geo- physics) to investigate the sub-surface down to 450 m depth (for locations of the boreholes see Fig. 6). In Fig. 7, simplified borehole lithologies are shown: Four boreholes have been drilled close to the church. While the main borehole Ky1/2014, (TLUG. LIAG), with ca. 457 m depth, crossed the Southern Kyffhäuser Fault, the two supple- mentary boreholes Ky2/2013, (TLUG) and Ky3/2013, (TLUG) only reached around 100 m in depth. The fourth borehole Ky1/1993, (TLUG), also around 100 m in depth, was drilled earlier. All four lithologies reveal around 3 m anthropogenic infill, then around 3-7 m Quarter- nary deposits. Below, the Permian Zechstein sequences start, with the Staßfurt-Anhydrite (z2/zS) from 7-74 m, the Staßfurt-Carbonate (Stinkschiefer) between 74-80 m depth. The full Werra-Anhydrite (z1/zW) is only pres- ent in Ky1/2014, reaching from ca. 80-160 m, followed by the Werra-Carbonate (Zechsteinkalk) between ca. 160-163 m. The Permian Zechstein sequences are com- pleted by the lowest marine deposit, the Kupferschiefer between ca. 163-165 m, followed by the final conglomer- ate layer. Below, between ca. 167-348 m, the Carbonifer- ous Kyffh äuser rocks are present (mainly sandstones). In 348 m depth, the Southern Kyffhäuser Fault is traversed, dipping with 50-70 o . Below, both Staßfurt-(z2/zS) and Werra-(z1/zW) deposits are found again, but this time with salt still present. Two older boreholes in vicinity of the sinkhole Quellgrund have been lowered in 1857 and 1866 to pro- vide brine water for medical purposes (e.g., Schriel & Bülow 1926a). The borehole Saline1 revealed a similar lithology as the main borehole Ky1/2014, traversing the Staßfurt and Werra rocks, then the Carboniferous Kyff- häuser rocks, crossing the Southern Kyffhäuser Fault, and finally the salt in the Staßfurt sequence below. The bore- Fig. 6: Simplified street map and topography of Bad Franken- hausen. Leaning church tower (red square), boreholes (blue stars with name of borehole), and Southern Kyffhäuser Fault (SKF, black dashed line). The Quellgrund is an ancient collapse sinkhole. Measured surface deformation is shown as red contours (vertical component, in mm) and as black lines with yellow arrow heads (horizontal components, longest arrow 90~mm). Fig. 7: Borehole lithologies, for leg- end see Fig. 3, For the more recent boreholes Ky*, cavities (white) and cavity infill (brown) from the cored samples is shown in the grey inlays in the right part of the lithologies). ACTA CARSOLOGICA 48/2 – 2019 181 GEORG KAUFMANN & DOUCHKO ROMANOV hole Saline2, south of the Southern Kyffhäuser Fault, re- veals a completely different lithology: Thick Tertiary and Triassic sequences til 217 m depth, then a thick Permian sequence (mainly Staßfurt anhydrite) down to the final depth of 423 m. In the recent boreholes, voids and cavities are also logged, which allows us to discuss the degree of karsti- fication in Bad Frankenhausen (Fig. 7). Voids logged in the four boreholes are frequent, ranging from heights of 20-50 cm up to the meter scale, with the largest ones ex- tending to more than 4 m vertically. Often, these voids are accompanied by collapse areas with infill from both brecciated remnants of gypsum and anhydrite, and resid- ual marls from the layered anhydrite. Most of the voids mapped are located in the Staßfurt-Anhydrit (z2/zS), with one void (1.45 m height) identified in the underly- ing Stinkschiefer (Staßfurt-Carbonate). The size of some of the voids indicates larger cavities, which might have similar dimensions to the mapped Schlotten. A georadar cross-borehole survey between boreholes Ky1/2014 and Ky1/1993 confirms these large-scale voids (Miensopust et al. 2015; Wadas et al. 2017). Thus we speculate that below the northern part of Bad Frankenhausen, larger sub-surface voids can be present and their partial break- down can be responsible for the observed surface defor- mations. FROM KARST TO CHEMISTRY … We have described the widespread karstification of the entire region, both along the outcropping Permian Zech- stein sequences and in the sub-surface. A special type of cavities, the Schlotten, entrance-less voids often discov- ered by mining operations, characterise large-scale hypo- genic karstification of the anhydrite layers, with aggres- sive water arriving from the carbonate layers below. In the next section, we discuss the chemical evolution of a solution on its way through the limestone into the anhy- drite above, with the potential to enlarge fractures to the meter-size and thus create large sub-surface entrance- less voids. CHEMISTRY In this section, we discuss the chemical properties of the soluble rocks involved (limestone, anhydrite, gypsum, salt) and the time scales of dissolution for these rocks types. We also consider the conversion of anhydrite to gypsum. REACTIONS The chemical reactions in the aqueous systems of gyp- sum (CaSO 4  • 2 H 2 O), anhydrite (CaSO 4 ), and limestone (CaCO 3 ) can be described as (e.g., Duan & Li 2008; Li & Duan 2011): (1) The first reaction describes the dissociation of water, the second the solution of carbon dioxide (CO 2 ) in water, the next two ones the dissociation of carbonic acid, and the remaining equations the solution of gypsum, anhydrite and limestone, respectively. EQUILIBRIA A solution being in contact with a surface of soluble rock will dissolve the rock, until a chemical equilibrium be- tween the dissolved species is attained. For gypsum, an- hydrite, and limestone, we can characterise this equilib- rium by tracking the calcium concentration, c [mol/m 3 ] and the respective calcium equilibrium concentration c eq [mol/m 3 ]. For gypsum and anhydrite, we can easily derive the calcium concentration from mass-balance consider- ations. For limestone, the derivation of an analytical re- lation for the calcium equilibrium concentration is a bit more tedious and results in an approximate solution (e.g., Dreybrodt 1988), which, however, compares well and to excellent accuracy with numerical solutions derived for the full electro-neutrality condition, which we verified by comparing our c eq for limestone to results obtained from PHREEQC (Parkhurst & Appelo 2013). The above mentioned chemical equations then re- sult in calcium equilibrium concentrations c eq for the dif- ferent soluble rock types: ACTA CARSOLOGICA 48/2 – 2019 182 MODELLING SPELEOGENESIS IN SOLUBLE ROCKS: A CASE STUDY FROM THE PERMIAN ZECHSTEIN SEQUENCES EXPOSED ALONG THE SOUTHERN HARZ MOUNTAINS AND THE KYFFHÄUSER HILLS, GERMANY (2) In (2), the mass-balance coefficients are: K A (T,z) (Blount & Dickson 1973) the equilibrium constant for the dissolution of anhydrite, K G (T,z) (Blount & Dickson 1973) the equilibrium constant for the disso- lution of gypsum, K H (T,z) (Weiss 1974) the equilibrium constant for the dissolution of atmospheric carbon dioxide into water (Henry’s law constant), K 0 (T,z) (Wissbrun et al. 1954) the equilibrium constant for the reaction of water and carbon dioxide to carbonic acid, K 1 (T,z) and K 2 (T,z) (Millero 1979; Mehrbach et al. 1973) the equilibrium con- stants for the dissociation of carbonic acid into bicarbon- ate, carbonate, and hydrogen, K C (T,z) (Mucci 1983) the equilibrium constant for dissolved calcite, γ Ca2+ , γ SO42- and γ HCO32- the activity coefficients for calcium, sulfate and bi- carbonate, p CO2 [atm] the carbon-dioxide partial pressure, T [ o C] the temperature of the solution, and z [m] is depth, accounting for hydraulic pressure increase with depth. For limestone, (2) is valid for the open system, in which the solution is in contact with the atmosphere, and carbon dioxide will be replenished by further solution of CO 2 from the atmosphere. However, most fracture en- largement takes place under closed-system conditions, and here the CO 2 is consumed and thus decreases with dissolution. In this case, the carbon-dioxide partial pres- sure is (Dreybrodt 1988): (3) with p atm CO2 [atm] the initial carbon-dioxide partial pres- sure obtained in the atmosphere and the soil. In Fig. 8, the calcium equilibrium concentrations are shown for limestone, anhydrite, and gypsum, and for comparison also the sodium equilibrium concentration for the dissolution of rock salt (NaCl) is shown as a function of temperature and for different depth (representing hydraulic pressure) and in the case of limestone also for different carbon-dioxide pressure values. While the equilibrium concentrations for limestone are between 0.1-3 mol/m 3 , values for gypsum and an- hydrite are about one order of magnitude higher, with values of around 15 mol/m 3 for gypsum and 40 mol/m 3 for anhydrite for temperatures below 40 o C. For halite, the sodium concentration is even two orders of magni- tude higher. Note that for temperatures below 50 o C, the solubility of gypsum is lower than the one for anhydrite, an important observation, which we will discuss in more detail below. Also note that for limestone and anhydrite, the equilibrium concentration is a retro-grade function Fig. 8: Calcium equilibrium con- centration as a function of temper- ature for limestone, anhydrite, and gypsum. Additionally, the sodium equilibrium concentration for halite is shown. Curves are shown for three different hydrostatic pres- sures (with hydraulic head as vari- able), and in the case of limestone for three different carbon-dioxide pressure representations. ACTA CARSOLOGICA 48/2 – 2019 183 GEORG KAUFMANN & DOUCHKO ROMANOV of the temperature, while for gypsum we observe a pro- grade dependence on temperature for temperatures be- low 30 o C. The dependence on hydraulic pressure is always below one order of magnitude, and higher pressures in- crease the equilibrium concentrations. For limestone, the dependence on carbon-dioxide pressure is non-linear and very pronounced, and respon- sible for the well-known and investigated dissolution processes in that rock. FLUX RATES With the calcium equilibrium concentration discussed, we now know, how much of the soluble rock surface can be dissolved. If we, however, also want to know, how fast the dissolution process occurs, we need to introduce the flux rate F [mol/m 2 /s]. Calcium flux rates are known both from labora- tory experiments and numerical modelling. For lime- stone, flux rates have been measured experimentally (e.g., Plummer et al. 1978; Svensson & Dreybrodt 1992; Eisenlohr et al. 1999), as well as determined by numerical experiments (e.g., Buhmann & Dreybrodt 1985a, 1985b; Dreybrodt & Kaufmann 2007; Kaufmann et al. 2010). For gypsum and anhydrite, flux rates have also been mea- sured in the laboratory (e.g., James & Lupton 1978; Go- bran & Miyamoto 1985; Lebedev & Lekhov 1990; Jeschke et al. 2001; Jeschke 2002). All of the above mentioned flux rate experiments can be described by a flux rates  F [mol/m 2 /s] with a piece-wise function of the calcium concentration c with respect to the calcium equilibrium concentration c eq : (4) with k i [mol/m 2 /s] a rate coefficient, m i [-] a flux rate regime coefficient, c [mol/m 3 ] the actual calcium concentration, c eq [mol/m 3 ] the calcium equilibrium concentration, and n i [-] a power-law exponent (see Kaufmann et al. (2010) for more details, and Tab. 1 for parameter values). The counter  i marks different re- gimes of the flux rate function. The rate coefficient de- pends on both surface-controlled rates and diffusion- controlled rates: (5) with k i S [mol/m 2 /s] the surface-controlled rate coefficient, and with k i D =2Dc eq /d [mol/m 2 /s] the diffusion-controlled rate coefficient, D [m 2 /s] the diffusion coefficient, and d [m] the film thickness. The film thickness describes the thickness of the solution layer on the fracture wall, e.g., for a circular conduit it can be defined by the ratio of flu- id volume in the conduit V [m 3 ] to reactive surface area of the conduit A (m 2 ), d=V/A. Tab. 1: Parameter values for soluble rock chemistry. Limestone atomic mass m Rock [kg/mol] 0.100 density ρ Rock [kg/m 3 ] 2600 rate constant 1 k 1 [mol/m 2 /s] 4x10 -7 rate constant k 2 =k 1 (1-c s ) n1-n2 rate exponent 1 n 1 [-] 1.0 rate exponent 1 n 2 [-] 4.0 Switch 1 c eq c s [-] 0.90 Anhydrite atomic mass m Rock [kg/mol] 0.136 density ρ Rock [kg/m 3 ] 2900 rate constant 2 k 1 [mol/m 2 /s] 4x10 -5 rate exponent 2 n 1 [-] 5.4 Gypsum atomic mass m Rock [kg/mol] 0.172 density ρ Rock [kg/m 3 ] 2200 rate constant 3 k 1 [mol/m 2 /s] 1.1x10 -3 rate constant k 2 =k 1 (1-c s ) n1-n2 rate exponent 3 n 1 [-] 1.0 rate exponent 3 n 2 [-] 4.5 Switch 3 c eq c s [-] 0.95 1 Buhmann et al. (1985a,b) 2 Jeschke (2002) 3 Jeschke et al. (2001) ACTA CARSOLOGICA 48/2 – 2019 184 MODELLING SPELEOGENESIS IN SOLUBLE ROCKS: A CASE STUDY FROM THE PERMIAN ZECHSTEIN SEQUENCES EXPOSED ALONG THE SOUTHERN HARZ MOUNTAINS AND THE KYFFHÄUSER HILLS, GERMANY In Fig. 9, flux rates following (4) are shown for an- hydrite, gypsum, and limestone. For both limestone and gypsum, the flux rate function is linear for strong under- saturation with respect to calcium, but becomes non-lin- ear close the the calcium equilibrium concentration. For anhydrite, the entire flux rate below saturation is non- linear. The non-linearity, which results in much lower flux rates close to saturation, is often a result of impuri- ties in the rock (e.g., shales, marls), which accumulate on the rock surface, while the soluble parts of the rock are dissolved and transported away in solution (Svensson & Dreybrodt 1992; Eisenlohr et al. 1999; Buhmann & Drey- brodt 1985a,b). Beyond the saturation point, precipitation of the mineral will start. Here, experimental data are rare, and only precipitation with calcite and limestone is based on experimental data (Buhmann & Dreybrodt 1985b; Drey- brodt & Buhmann 1991). For precipitation rates of anhy- drite and gypsum, we simply assume a linear relation (see Kaufmann et al. 2017 for more details). This assumption rests on the fact, that during dissolution impurities ac- cumulate on the interface solution-rock, and are respon- sible for the non-linear rates close to equilibrium, while during precipitation the pure mineral species will depos- it, which we assume can be described by a linear rate law. TIME SCALES We now know what the amount of dissolved material is (c eq ) and how fast is will be dissolved (F). We combine these two informations to describe the temporal evolu- tion of the calcium concentration c in the water film with a differential equation: (6) with A s [m 2 ] the reactive surface area, and V [m 3 ] the fluid volume. The ratio A s /V is known as specific surface, and we can express this in terms of film thickness d by defin- ing a circular conduit with surface area A s =2 p  r l and vol- ume V= p  r 2  l, with r [m] the conduit radius (d=2r) and l [m] the conduit length. We thus find A s /V~1/r. Similar relations hold for rectangular conduits. Note that (6) de- scribes a stagnant film, transport and advection will be considered in the next section. With both c eq and F known for anhydrite, gypsum, and limestone, the temporal evolution of the modelled calcium concentration c is shown in Fig. 10. All solutions start with a calcium concentration of zero, c in =0. For the low film thickness value of d=0.001  m (solid lines), a water film on limestone reaches equilibrium after about 10 h, while the same solution on a gypsum surface would be saturated after only 1 h. For anhydrite, the non-linear flux rate is responsible for reaching the saturation point even later, around 50 h in this case. Note the curve for halite, which we plotted for comparison. The sodium equilibrium concentration in this case is reached within minutes! If we increase the film thickness by an order of magnitude (d=0.010 m, dotted lines), the times needed to reach the equilibrium concentrations for the soluble Fig. 9: Flux rates as a function of calcium concentration normal- ized to calcium equilibrium con- centration. The displayed flux rates are calculated for a tem- perature of T=10 o C, a hydraulic head of z=0 m, a film thickness of d=0.00001 m, and a carbon-diox- ide pressure of p CO2 =0.05 atm for calcite. Note the different scaling constants for the flux rates for the different rock types. ACTA CARSOLOGICA 48/2 – 2019 185 GEORG KAUFMANN & DOUCHKO ROMANOV rock types shown increase by about one order of mag- nitude. The reason for the increase in time with increas- ing fracture width is the boundary layer, which for this larger fracture is more pronounced and the diffusion of dissolved species out of this boundary layer takes longer (e.g., Buhmann & Dreybrodt 1985a, 1985b). Thus, in- creasing the film thickness d changes the rate coefficient (5). ANHYDRITE-GYPSUM CONVERSION We now finally discuss a water film in a fracture in anhy- drite, with solution undersaturated with respect to cal- cium passing through it. As discussed above, the calcium concentration c increases in the fracture, but now we al- low gypsum to precipitate, once the calcium equilibrium concentration for gypsum, c eq gypsum ≈15 mol/m 3 , is reached (Fig. 11). Note that the water film is still stagnant. Fig. 10: Calcium concentration as a function of time. Shown are con- centrations for anhydrite, gypsum, and limestone approaching their respective calcium equilibrium concentrations (dashed lines), for a temperature of T=10 o C, a hydrau- lic head of z=0 m, a film thickness of d=0.001  m (solid lines) and a film thickness of d=0.01 m (dotted lines), and a carbon-dioxide pres- sure of p CO2 =0.05  atm for calcite. Additionally, for halite the sodium concentration and the respective sodium equilibrium concentration is also shown (note that these quan- tities are scaled by a factor 1/200)! Fig. 11: Calcium concentration as a function of time in the case of anhy- drite-gypsum conversion. Once the calcium concentration of the solu- tion passes the calcium equilibrium concentration of gypsum, gypsum starts to precipitate. Also shown is the reduction in fracture width d as a result of the gypsum precipitation for two initial fracture width (see legend). ACTA CARSOLOGICA 48/2 – 2019 186 MODELLING SPELEOGENESIS IN SOLUBLE ROCKS: A CASE STUDY FROM THE PERMIAN ZECHSTEIN SEQUENCES EXPOSED ALONG THE SOUTHERN HARZ MOUNTAINS AND THE KYFFHÄUSER HILLS, GERMANY We follow a strategy first described by Kontrec et al. (2002) and later extended by Serafeimidis and Anagnos- tou (2012). We solve the coupled reaction equations for anhydrite and gypsum concentration and fracture width: (7) In (7), r [m] is the fracture half-width, and m anh,gyp [kg/mol] the atomic mass and ρ anh,gyp [kg/m 3 ] the density of anhydrite or gypsum, with anhydrite the dissolving and gypsum precipitating mineral. Note that the terms for the gypsum flux rate, F gyp , are only active, when the calcium concentration c in the solution exceeds the cal- cium equilibrium concentration for gypsum, c eq gypsum , be- cause at the start of the experiment no gypsum is present. The concentration curve for anhydrite dissolution starts as before (green line), increasing from zero and heading towards the calcium equilibrium concentration of anhydrite, c eq anhydrite »45 mol/m 3 . However, soon after pass- ing c eq gypsum ≈15  mol/m 3 (blue dashed line), gypsum starts to precipitate and the concentration curve approaches a plateau at c≈17-18 mol/m 3 . With time, the concentration slowly drops again. The reason for this drop is that gyp- sum precipitates faster than anhydrite dissolves, and the fracture width therefore slowly reduces. This can be seen for the two black lines with different initial fracture widths (d=0.001  m, solid black line, d=0.0005  m, dotted black line): Both fractures initially grow a bit in size (hardly vis- ible in the plot), but once the calcium concentration ex- ceeds the calcium equilibrium concentration for gypsum (c>c eq gypsum ), the fracture width decreases. This closure of the fracture is forced by the larger volume the gypsum oc- cupies, when compared to anhydrite. After about 100 h (roughly four days) for the frac- ture with d=0.0005  m initial width and about 1000 h (roughly a month) for the fracture with d=0.001 m initial width, the clogging accelerates and after around 1 year the fracture is clogged. Only fractures with larger initial width d>0.01  m remain open, as from this size on the amount of gypsum precipitating is no longer enough to fill the fracture completely. FROM CHEMISTRY TO EVOLUTION … We have described the dissolution (and precipitation) of soluble rocks for the three important rock types lime- stone, gypsum, and anhydrite present in the Permian Zechstein sequences in the Southern Harz and the Kyff- häuser Hills. For our geological setup, it is important that meteoric water is transported deeper into the sub-surface through the well karstified carbonates, often saturated with respect to this carbonates, but still able to dissolve anhydrite, which has a much higher calcium equilibrium concentration. As we are most likely still in a flow re- gime with temperatures below 60 o C, during the dissolu- tion of anhydrite the calcium equilibrium concentration for gypsum is reached, and gypsum starts to precipitate. Depending on the initial fracture size distribution in the anhydrite, smaller fractures are quickly clogged and then seal off the anhydrite (as aquitard). Only a thin layer of anhydrite along the carbonate-anhydrite boundary is dissolved, producing the Schlottenlehm, a residual clay material often found along this interface. However, in areas characterised by larger initial frac- tures, e.g., faults, the precipitated gypsum cannot clog the entire fracture, and the anhydrite dissolution proceeds deeper into the anhydrite rocks, creating the initial voids needed for the evolution of the Schlotten. We speculate that the observed surface deformation with its visible landmark the tilted church tower of the Oberkirche in Bad Frankenhausen is also a result of solution processes within the Permian Zechstein rocks beneath the city. With this knowledge of dissolution and precipita- tion in the soluble Permian Zechstein rocks, we proceed to the next section, discussing a three-dimensional mod- el for the evolution of flow and transport in the system. We aim to describe a simplified system with hypogene water supply through limestone into anhydrite, and the creation of local cavities, the Schlotten. EVOLUTION In this section, we describe the numerical tool SOCKS3D (SOluble roCKS), which is used to describe the evolu- tion of flow and permeability in an aquifer, which is in parts composed of soluble rocks. Our aim is to describe a situation as depicted in Fig. 12: The hypothetical model describes a simplified cross section similar to the Bad Frankenhausen region, with a single sequence of Per- mian Zechstein rocks (limestone, anhydrite, sandstone), offset by a fault. On the left side of the fault, the soluble rocks are in greater depth, and here salt is present on top of the anhydrite. ACTA CARSOLOGICA 48/2 – 2019 187 GEORG KAUFMANN & DOUCHKO ROMANOV MATHEMATICAL FRAMEWORK We use the numerical program package SOCKS3D, which is a new fork of our widely used KARSTAQUI- FER package (see Kaufmann et al. 2010; Hiller et al. 2011; Kaufmann et al., 2012; Hiller et al. 2014, for more details), including the description of multiple rock types. The KARSTAQUIFER package has been used to simulate the speleogenetic evolution of an aquifer comprised of soluble and/or insoluble rock. Here, we will provide only a short model description, the reader interested in more details is referred to the references listed. The model is based on a two-dimensional digital to- pography (nx x ny grid points), and then is extended into depth with nz-1-layers, resulting in nx x ny x nz nodal points. On the discretisation level, rock matrix is iden- tified by 3D parallelepipetal elements, rock fractures by 1D linear bar elements, both represented with finite ele- ments. The SOCKS3D packages solves the following set of equations in a three-dimensional modelling domain describing a fractured, porous aquifer in soluble rocks under steady-state or transient flow conditions and a re- alistic topography: (8) The first equation is a conservation equation for flow, with h [m] the hydraulic head, [m/s] the Darcy veloc- ity, S s [1/m] the specific storage, [1/m] the nabla op- erator, and t [s] time. The second equation is the trans- port equation for the calcium concentration c [mol/m 3 ], with P=πd [m] and A=πd 2 /4 [m 2 ] the perimeter and the cross section of a circular fracture (not to confuse with A s , the specific surface area), d [m] the fracture diameter, F [mol/m 2 /s] the calcium flux rate, and c [mol/m 3 ] the calcium concentration. The third equation describes the reactive enlargement of a fracture with m rock [kg/mol] the atomic mass and ρ rock [kg/m 3 ] the density of the soluble rock. Darcy velocity is given by (9) with K [m/s] the hydraulic conductivity. For fracture elements, the hydraulic conductivity is given as (10) with ρ [kg/m 3 ] the density of the fluid, g [m/s 2 ] the gravi- tational acceleration, η [Pa s] the viscosity of the fluid, f [-] the friction factor, and q i [m/s] the Darcy flow rate component in flow direction from the previous itera- tion step. Note that flow in the presence of turbulence becomes non-linear. For matrix elements, the hydraulic conductivity is derived from the Kozeny-Carman-rela- tion (11) with F  [-] the porosity of the rock matrix, scaled between 0 (no pores) and 1 (only pores), l * a tortuosity factor ac- counting for the complex geometry of the flow through the pores, and d’ defined as typical grain size. MODEL SETUP The three-dimensional model extends 100 m in north- south and 200 m in west-east direction, and 20 m verti- cally, using 51x101x11 nodal points (Fig.13). With this setup, we simulate a section of the Zechstein deposits of the southern Harz Mountains, simplified with one depo- sitional sequence. The model has a 2  m resolution, with around 160,000 fracture and 50,000 matrix elements. While the hydraulic properties of the matrix elements are fixed to a hydraulic conductivity of K m ≈10 -7   m/s and a specific storativity of S m =5x10 -6  1/m, hydraulic conductivities of Fig. 12: Hypothetical model describing the conditions along a se- quence of Zechstein rocks (color coded) offset by faulting (black line). Water enters the system either as meteoric water via precipi- tation (along the exposed soluble rocks), or as hydrothermal water ascending through the fault zone. ACTA CARSOLOGICA 48/2 – 2019 188 MODELLING SPELEOGENESIS IN SOLUBLE ROCKS: A CASE STUDY FROM THE PERMIAN ZECHSTEIN SEQUENCES EXPOSED ALONG THE SOUTHERN HARZ MOUNTAINS AND THE KYFFHÄUSER HILLS, GERMANY the fracture elements change with time due to dissolution and precipitation of the mineral species. Specific storativ- ity of the fracture is low, S f =10 -8  1/m. The model consists of up to three different rock types, from top to bottom sandstone, anhydrite, and limestone. Surface temperatures are set to T=10 o C. As boundary conditions, we define a fixed head boundary condition along the left side (h=10  m, blue dots), defining a regional base level of the aquifer sys- tem. As input, a fixed inflow condition along the bottom part of the model along the right side is defined, mim- icking a regional groundwater-flow component into the limestone layer (q=20000 m 3 /yr, T=10 o C, p CO2 =0.02 atm, c eq =1.01 mol/m 3 , c in =0.90 mol/m 3 , red dots). Thus the in- flow occurs into the limestone layer, while the local base level is the gypsum/anhydrite formation along the right side. This conditions mimic hypogene conditions for the evolution of cavities. The resulting initial groundwater-flow pattern can be seen in Fig. 13: The head contours reveal flow from right to left, with only small local perturbations resulting from the initially heterogeneous fracture conductivities visible in the distorted head contours, which result from our assignment of initial fracture diameter values using a log-normal distribution, assigning an average value d ave  [m] and a standard deviation d σ [m]. MODEL RUNS With the geometrical, hydrological, and chemical prop- erties defined, we discuss a set of model runs for the evolving aquifer, which are all run for 20,000 years. Our strategy is as follows: We start with a simple one-material limestone aquifer below a sandstone domain to realise how the carbonates of the Zechstein sequences become karstified, both in time and geometrically. We then add a gypsum layer between the limestone and the sandstone to understand the impact of dissolution of different sol- uble rock types, resulting in different time scales of the evolution of secondary permeability in both soluble rock sequences. We then swap the gypsum layer with an an- hydrite layer to describe the simplified setup of the Zech- stein sequence in the southern Harz region. Our aim is to understand the clogging of the anhydrite sequence due to the precipitation of secondary gypsum. As we will see, we efficiently block the anhydrite, and therefore we finally need to introduce conditions, which can lead to Fig. 13: Setup of numerical evo- lution model. Top: Schematical cross section. Bottom: Topography (transparent), limestone (red), an- hydrite (grey), sandstone (blue). Inflow boundary (right, red dots), and base level (left, blue dots). Groundwater-flow head contours (colour-coded). The block in the background marks the different rock layers (for colour code see Fig. 13 top). ACTA CARSOLOGICA 48/2 – 2019 189 GEORG KAUFMANN & DOUCHKO ROMANOV localised enlarged voids in the anhydrite, resembling the Schlotten-type cavities. LIMESTONE AQUIFER The first model is made up of a 6 m thick limestone layer and a 14 m thick sandstone rock above. With this model, we introduce a standard model of karst aquifer evolution with its typical time scale of evolution. The sandstone has fairly low permeability, and re- charge is only via a regional flow component from the right into the limestone layer at depth, mimicking the regional flow component. The initial fracture diameter distribution for both the sandstone and the limestone is based on a log-normal distribution with d ave =0.1 mm and d s =0.05 mm. With this large standard deviation, we expect several competing preferential flow paths, as the distribution of initial frac- ture diameter values becomes wide enough to provide less direct flow path options. Water entering the limestone layer from the right side is slightly aggressive, with c in =0.90  mol/m 3 it has already attained 90\% of the equilibrium values of c eq =1.01  mol/m 3 for this temperature / carbon-dioxide pressure setting for limestone. The evolution of this mod- el is shown in Fig. 14 for two snapshots in time. 5,000 years into the evolution, several clusters of enlarged fractures have started developing from the right side towards the base level at the left side. En- larged fractures have already diameter values in the cm- range. The different clusters compete for flow, and the bottommost cluster has already attained an advantage, when compared to the other clusters, which can be seen from the size of the cluster and the head distribution in front of it. Fig. 14: Limestone aquifer: Two snapshots in time, with fracture di- ameter and hydraulic head colour- coded. The block in the background marks the different rock layers (for colour code see Fig. 13 top). ACTA CARSOLOGICA 48/2 – 2019 190 MODELLING SPELEOGENESIS IN SOLUBLE ROCKS: A CASE STUDY FROM THE PERMIAN ZECHSTEIN SEQUENCES EXPOSED ALONG THE SOUTHERN HARZ MOUNTAINS AND THE KYFFHÄUSER HILLS, GERMANY 20,000  years into the evolution, the bottommost cluster has succeeded in capturing the majority of flow, and a preferential flow path connecting the right and the left side is established, with fracture diameter values in the 10 cm-range. Hydraulic pressure drops signifi- cantly, and the further evolution is deferred due to the lower flow velocities in this stage. This situation, called catchment-controlled evolution in the karst aquifer litera- ture, limits flow through the evolving aquifer, as inflow is controlled by the amount of water coming from a limited catchment. LIMESTONE-GYPSUM AQUIFER We proceed with the second model, which is made up of a 6 m thick limestone layer, a 6 m thick gypsum layer, and a 8 m thick sandstone rock above. With this second model, we introduce a model of karst aquifer evolution made up of two adjacent layers of different soluble rock. The different properties of the limestone and gypsum rocks will cause a completely different evolution, as we will see below. The model is recharged as above, via the right side of the limestone part, with water close to the calcium sat- uration for limestone under the chosen climatic condi- tions (c in limestone =0.90 mol/m 3 , c eq limestone =1.01 mol/m 3 ). Thus the incoming water can dissolve limestone (in the non- linear regime of the limestone flux rate), but has much more potential to dissolve the overlying gypsum, which with the chosen surface temperature has a calcium equi- librium concentration of c eq gypsum =14.57 mol/m 3 . The evo- lution of this model is shown in Fig. 15 for two snapshots in time. After already 300 years of evolution, the enlarge- ment of fractures in the gypsum aquifer is substantial, Fig. 15: Limestone-Gypsum aqui- fer: Two snapshots in time, with fracture diameter and hydraulic head colour-coded. The block in the background marks the differ- ent rock layers (for colour code see Fig. 13 top). ACTA CARSOLOGICA 48/2 – 2019 191 GEORG KAUFMANN & DOUCHKO ROMANOV with clusters enlarged to the 10 cm size. After 20,000 years of evolution, the entire gypsum sequence is basically dis- solved. Hardly visible from the current perspective view and therefore only described is the fact, that in the lime- stone sequence almost no evolution of secondary poros- ity occurs. The reason for this behaviour is that water in- jected into the limestone along the right size of the model domain finds an easier flow path through the gypsum se- quence, which through its fast dissolution becomes much more permeable than the limestone and attracts flow. Thus, the dissolution kinetics (gypsum versus limestone) entirely controls the evolution of this combined system of two soluble rocks, with gypsum having a clear advantage due to its much higher solubility. LIMESTONE-ANHYDRITE AQUIFER We, however, know that in the southern Harz Mountains the limestone is capped by anhydrite. While both gypsum and anhydrite are made up of the same material, consid- ering the evolution of a combined limestone-anhydrite system is completely different, as we will see below in the snapshot shown in Fig. 16. In this model, we keep the hydraulic boundary con- ditions as above, but in the anhydrite layer we allow for gypsum precipitation, once the calcium concentration passed the calcium equilibrium concentration for gyp- sum, c eq gypsum ≈15.34 mol/m 3 . As we have already discussed in section 4, this combined anhydrite dissolution-gyp- sum precipitation can clog small fractures on time scales below a year. Thus it is not surprising to find an evolved aquifer, which after 20,000 years of evolution time re- sembles more or less the situation of the single limestone aquifer: A couple of enlarged clusters of fractures in the limestone part of the model, with the bottommost cluster on its way to win the competition. However, no evolution of secondary porosity in the anhydrite sequence, as here the water entering from the limestone layer below starts dissolving the anhydrite, but quickly passes the gypsum equilibrium concentration, resulting in gypsum pre- cipitation. As the precipitated gypsum occupies a larger volume as the dissolved anhydrite, the small fractures initially present in the anhydrite clog and the anhydrite layer becomes almost impermeable. With this model we have shown, how effective the anhydrite-gypsum conversion is in reducing the perme- ability of an anhydrite layer, resulting in properties more characteristic for an aquitard for the anhydrite, which is a common observation in nature. LIMESTONE-ANHYDRITE AQUIFER WITH SCHLOTTEN-TYPE CAVITIES Now we have to come back to our initial question, how do Schlotten-like cavities evolve in a combined lime- stone-anhydrite sequence as observed in the southern Harz Mountains. As we are able to clog the anhydrite recharged with an aggressive water source from below (hypogene setting), we need some conditions to allow for the evolution of secondary porosity locally. We have sev- eral options available, which can be discussed with our knowledge from section 4. Two options applicable to the southern Harz Mountains and the Kyffhäuser setting will be discussed: • C he mi c a l c o n t r o l: Coming back to Fig. 8, we have based our discussion of the effective clogging of the anhydrite layer on the different calcium equilibrium concentrations of gypsum and anhydrite, with around c eq gypsum ≈15 mol/m 3 and c eq anhydrite »45 mol/m 3 for a tem- perature of 10 o C. This large difference in c eq and the volume increase between the dissolved anhydrite and Fig. 16: Limestone-Anhydrite aq- uifer: One snapshot in time, with fracture diameter and hydraulic head colour-coded. The block in the background marks the differ- ent rock layers (for colour code see Fig. 13 top). ACTA CARSOLOGICA 48/2 – 2019 192 MODELLING SPELEOGENESIS IN SOLUBLE ROCKS: A CASE STUDY FROM THE PERMIAN ZECHSTEIN SEQUENCES EXPOSED ALONG THE SOUTHERN HARZ MOUNTAINS AND THE KYFFHÄUSER HILLS, GERMANY the gypsum precipitate causes the clogging of the layer. However, when we increase the temperature of the so- lution to 30-40 o C, this difference becomes very small, and the gypsum precipitation less efficient. We argue that ascending thermal water through a fault zone (as the southern Kyffhäuser Fault in Bad Frankenhausen) can cause this reduction in gypsum precipitation, thus leaving voids from enlarged anhydrite fissures open. • S t r u c t ur a l c o n t r o l: Another possible mechanism to create localised voids in the anhydrite can be derived from an observation in Fig. 11: The time scale for clog- ging an anhydrite fracture with precipitating gypsum depends on the initial width of the anhydrite fracture, larger initial fractures will clog later. This effect can be substantial, as small increases in initial fracture width already result in orders-of-magnitude delay in clog- ging the fracture. We argue that a fault through the Zechstein sequence, providing a larger initial perme- ability of the rock sequences, can allow fractures in the anhydrite to remain open and to enlarge to larger voids. We will show that both hypotheses result in local- ised enlarged voids in the anhydrite layer in our simpli- fied model. The model shown in Fig. 17 is based on the limestone-anhydrite-sandstone setup of Fig. 16, but with two vertical fault zones (shown as vertical semi-transpar- ent planes). The fault on the left side close to the base level has an initial fracture diameter of 0.001 m, with is one order of magnitude larger than the average initial fracture di- ameter of the other fractures (assigned with log-normal distribution). The fault on the right side of the model domain has an initial fracture diameter of 0.005 m, thus it is hydraulically more effective, when compared to the Fig. 17: Limestone-Anhydrite aqui- fer with Schlotten-type cavities: Two snapshot in times, with frac- ture diameter and hydraulic head colour-coded. The two faults are shown as vertical semi-transparent planes. In the bottom figure, only fractures enlarged to the 50 cm size and above are shown. The block in the background marks the differ- ent rock layers (for colour code see Fig. 13 top). ACTA CARSOLOGICA 48/2 – 2019 193 GEORG KAUFMANN & DOUCHKO ROMANOV other fault. While the fault on the right side is used as structural-control element, along the fault on the left side we inject water from below with a temperature of 40 o C, p CO2 =0.02  atm, and a calcium concentration of c eq =1.01  mol/m 3 . Thus, this thermal water component is saturated with calcium with respect to the limestone chemical-control element. Pointing back to our Fig. 12, the left fault mimics groundwater ascending through the southern Kyffhäuser Fault in Bad Frankenhausen as the origin of the enlarged voids observed in the boreholes, while the right fault depicts a situation similar to the Elis- abethschächter Schlotte with its relation to the collapse sinkhole on the surface. 10,000 years into the evolution, we observe a pref- erential enlargement of fractures in the limestone part of the rock sequence, but the pattern is different when compared to our other runs. The reason is the right fault, which attracts flow because of its larger initial fracture diameter values. The fault acts as local base level, from which new preferential pathways in the limestone start to evolve towards the base level on the left side. In the an- hydrite part of the two faults, several anhydrite fractures remain open and enlarge width time to the cm-to-meter range. In the second snapshot taken 20,000 years into the evolution, we just visualize fractures enlarged to at least 50 cm to clarify the substantial local growth in the anhy- drite layer. In the anhydrite part of the right fault, several an- hydrite fractures remain open and become enlarged to the meter size, generating a local system of intercon- nected enlarged voids, similar to the Elisabethschächter Schlotte! The evolution around the left fault, in which we inject thermal water, is different. Here, several isolated clusters of enlarged voids grow on an downstream of the fault, which quickly reach sizes in the 10 cm to meter range. Such voids are possibly observed in the boreholes lowered around the church Oberkirche in Bad Franken- hausen! We note, however, that here also water from the neighbouring sandstone can provide additional aggres- sive solution, as in the case of the Numburg Cave. CONCLUSIONS In this final section, we review our models and summa- rize our key results in view of the questions asked in the introduction. We have asked three key questions in the introduction, which we now can answer: 1. Can we explain the evolution of the Schlotten-type cavi- ties with a large-scale karst evolution model, account- ing for dissolution in a combined carbonate-sulfate setting? Yes, we can. Our numerical modelling results suggest that in a setting where water seeps into the Permian Zechstein sequences through the carbonate layers (Zechsteinkalk, Stink\-schiefer), it dissolves the anhydrite from below. However, the dissolution of the anhydrite is mainly occurring along the inter- face carbonate-anhydrite. Fissures in the anhydrite are usually clogged by gypsum precipitation. Only along prominent fault zones, flow is diverted and focussed into these hydraulically more permeable zones and anhydrite dissolution outpaces the gypsum precipita- tion. Thus, along these prominent fault zones, Schlot- ten-type cavities can form. 2. Can we explain the surface deformation observed in Bad Frankenhausen with numerical models of large- scale dissolution in the underlying sulfate layers? Yes, we can, at least indirectly through the evolution of large isolated cavities in near-surface layers along prominent fault zones, which develop in the anhydrite and due to their sizes (often in the meter range, as our models suggest), these voids are prone to collapse. 3. Can we estimate the role of aggressive mineral-rich brine water arriving from below through fault systems (hypogene setting)? We can also speculate that rising thermal water (as in our setup of one of the prominent faults) changes the chemical equilibrium conditions and thus enhances the dissolution of anhydrite, result- ing in even faster development of large sub-surface voids. ACTA CARSOLOGICA 48/2 – 2019 194 MODELLING SPELEOGENESIS IN SOLUBLE ROCKS: A CASE STUDY FROM THE PERMIAN ZECHSTEIN SEQUENCES EXPOSED ALONG THE SOUTHERN HARZ MOUNTAINS AND THE KYFFHÄUSER HILLS, GERMANY ACKNOWLEDGMENT We would like to thank Wolfgang Dreybrodt for his very constructive review. We acknowledge funding from the BMBF within the SIMULTAN project under research grant 03G0843G (GK) and from the DFG under research grant KA1723/6-2 (DR). Figures were prepared using GMT software (Wessel & Smith 1998, Wessel et al. 2013). Borehole logs are kindly provided by TLUG (Thüringer Landesanstalt für Umwelt und Geologie). REFERENCES Adler, A., 2013: Die känozoische Entwicklung des Gewäss- ersystems der Barbarossahöhle (Kyffäuser): Hydro- geologische und paläokologische Untersuchungen.- Master’s thesis, Martin-Luther-Universität Halle- Wittenberg, Halle, pp. 242. Adler, A. & D. Mertmann, 2012: Fossilien im Höhlen- lehm der Barbarossahöhle, Kyffhäuser.- Mitt. Verb. Dt. Höhlen- und Karstforscher, 58 (1), 6-16. Bauer, S. & A. Elste, 2016: Karst und Bergbau im Südöstlichen Harzvorland, Sachsen-Anhalt. In: Grubenarchäologische Gesellschaft e.V. (eds.) 19. Internationaler Bergbau & Montanhistorik Work- shop Mansfeld-Südharz 2016. Clausthal-Zellerfeld, pp. 95-107. Blount, C. & F. Dickson, 1973: Gypsum-anhydrite equi- libria in systems CaSO 4 -H 2 O and CaSO 4 -NaCl- H 2 O.- Am. Mineralogist, 58, 323-331. Brust, M., 2005: Die Erschließung und Entwicklung der Barbarossahöhle als Schauhöhle.- In: Kupetz, M., Brust, M. (eds.) Karst und Altbergbau am Kyffhäus- er, Salz, Kupfer, Gips, Alabaster. - Tagungspublika- tion zum 17. Treffen des Arbeitskreises Bergbau- folgelandschaften vom 08.-09.04.2005 im Geopark Barbarossahöhle. Hannover, 225, pp. 18-23. Brust, M., 2010: Neuer Erdfall nahe der Äbtissinnengrube bei Bad Frankenhausen 2009 - Kurzbericht und Analyse der Medienberichterstattung.- Mitt. Ver- band dt. Höhlen- und Karstforscher, 56(4), 113-115. Brust, M. & J. Graf, 2016: Die Schlotte auf dem Segen- Gottes-Stollen und die Elisabethschächter Schlotte.- Grubenarchäologische Gesellschaft [Online] Avail- able from: http://www.untertage.com/downloads/ Brust_et_Graf_2016_Die_Schlotte_auf_dem_Se- gen-Gottes-Stollen_etc.pdf. Brust, M. & G. Nash, 2013: Mine caves on the south-east- ern flank of the Harz Mountains (Saxony-Anhalt, Germany).- In: De Waele, J., Forti, P., Naseddu, A. (eds.) Mine Caves: Grotte di Miniera. Vol. Ser. II, Vol. XXVIII. pp. 71-72. Buhmann, D. & W. Dreybrodt, 1985a: The kinetics of calcite dissolution and precipitation in geologi- cally relevant situations of karst areas. 1. Open system.- Chem. Geol., 48, 189-211, https://doi. org/10.1016/0009-2541(85)90046-4. Buhmann, D. & W. Dreybrodt, 1985b: The kinetics of calcite dissolution and precipitation in geologi- cally relevant situations of karst areas. 2. Closed system.- Chem. Geol., 53, 109-124, https://doi. org/10.1016/0009-2541(85)90024-5. De Waele, J., Forti, P . & A. Naseddu, 2013: Introduction.- In: De Waele, J., Forti, P., Naseddu, A. (eds.) Mine Caves: Grotte di Miniera. Vol. Serie II, Vol. XXVIII. pp. 11-16. Dreybrodt, W ., 1988: Processes in Karst Systems.- Spring- er, pp. 288, Berlin, https://doi.org/10.1007/978-3- 642-83352-6. Dreybrodt, W. & D.  Buhmann, 1991: A mass trans- fer model for dissolution and precipitation of cal- cite from solutions in turbulent motion.- Chem. Geol., 90, 107-122, https://doi.org/10.1016/0009- 2541(91)90037-R. Dreybrodt, W. & G.  Kaufmann, 2007: Physics and chemistry of dissolution on subaerially exposed soluble rocks by flowing water films.- Acta Car- sologica, 36 (3), 357-367, https://doi.org/10.3986/ ac.v36i3.169. Duan, Z. & D. Li, 2008: An improved model for the calculation of CO 2 solubility in aqueous solutions containing Na + , K + , Ca 2+ , Mg 2+ , Cl - , and SO 4 2- .- Geo- chem. Cosmochem. Acta, 72, 5128-5145, https:// doi.org/10.1016/j.marchem.2005.09.001. Eisenlohr, L., Meteva, K., Gabrovšek, F. & W.  Drey- brodt, 1999: The inhibiting action of intrinsic impurities in natural calcium carbonate miner- als to their dissolution kinetics in aqueous H 2 O- CO 2 solutions.- Geochim. Cosmochem. Acta, 63 (6), 989-1001, https://doi.org/10.1016/S0016- 7037(98)00301-9. Gobran, G. & S. Miyamoto, 1985: Dissolution rates of gypsum in aqueous salt solutions.- Soil Science, 140 (2), 89-93, https://doi.org/10.1097/00010694- 198508000-00002. ACTA CARSOLOGICA 48/2 – 2019 195 GEORG KAUFMANN & DOUCHKO ROMANOV Hiller, T., Kaufmann, G. & D. Romanov, 2011: Karstifica- tion beneath dam sites: From conceptual models to realistic scenarios.- J. Hydrol., 398, 202-211, https:// doi.org/10.1016/j.jhydrol.2010.12.014. Hiller, T., Romanov, D., Gabrovšek, F. & G. Kaufmann, 2014: The creation of collapse dolines: a 3D model- ling approach.- Acta Carsologica, 43 (3-4), 241-255, https://doi.org/10.3986/ac.v43i2.832. Hohl, R., 1985: Die Entwicklungsgeschichte der Erde.- Verlag Werner Dausien, Hanau/Main, 6 th edition, pp. 703. James, A. & A. Lupton, 1978: Gypsum and anhydrite in foundations of hydraulic structures. Geotech- nique, 28 (3), 249-272, https://doi.org/10.1680/ geot.1978.28.3.249. Jeschke, A., Vosbeck, K. & W. Dreybrodt, 2001: Surface controlled dissolution rates of gypsum in aqueous solutions exhibit nonlinear dissolution rates.- Geo- chem. Cosmochem. Acta, 65 (1), 27-34, https://doi. org/10.1016/S0016-7037(00)00510-X. Jeschke, A. A., 2002: Lösungskinetik von Gips und An- hydrit.- Ph.D. thesis, Universität Bremen, Bremen, Germany. Kaufmann, G., Gabrovšek, F. & D. Romanov, 2017: Frac- ture evolution in soluble rocks: From single-ma- terial fractures towards multi-material fractures.- Acta Carsologica, 46/2-3, 199-216, https://doi. org/10.3986/ac.v46i2-3.5135. Kaufmann, G., Romanov, D. & W.  Dreybrodt, 2012: Modeling of karst aquifers.- In: White, W., Culver, D. & Pipan, T. (eds.) Encyclopedia of Caves. Aca- demic Press, Elsevier, 3 rd edition, pp. 1025. Kaufmann, G., Romanov, D. & T. Hiller, 2010: Model- ling three-dimensional karst aquifer evolution using different matrix-flow components.- J. Hy- drol., 388, 241-250, https://doi.org/10.1016/j.jhy- drol.2010.05.001 Kempe, S., 2008: Gipskarst - ein Überblick.- In: Brust, M. K., Kupetz, M. & S. Schmiedel (eds.) Gips- und Anhydritkarst in der Mansfelder Mulde – Die Wim- melburger Schlotten, 23. Treffen Arbeitskreis Berg- baufolgen, Deutsche Gesellschaft für Geowissen- schaften. 235, pp. 30-41. Kempe, S., Bauer, I. & S. Glaser, 2017: Hypogene caves in Germany.- In: Klimchouk, A., N. Palmer, A., De Waele, J., Auler, A.S. & P. Audra (eds.) Hypogene Karst Regions and Caves of the World. Springer, pp. 329-347. Kersten, T., Kobe, M., Gabriel, G., Timmen, L., Schön, S. & D. Vogel, 2017: Geodetic monitoring of sub- rosion-induced subsidence processes in urban ar- eas.- J. Appl. Geodesy, 11  (1), 21-29, https://doi. org/10.1515/jag-2016-0029. Klimchouk, A., 2011: Hypogene Speleogenesis: Hydro- geological and Morphogenetic Perspective.- NCKRI Special Paper No. 1. Kontrec, J., Kralj, D. & L. Brecevic, 2002: Transforma- tion of anhydrous calcium sulphate into calcium sulphate dihydrate in aqueous solutions.- Jour- nal of Crystal Growth, 240, 203-211, https://doi. org/10.1016/S0022-0248(02)00858-8. Kupetz, M., 2005: Gipskarst am Kyffhäuser und die Gen- ese der Barbarossahöhle.- In: Kupetz, M., Brust, M. (eds.) Karst und Altbergbau am Kyffhäuser, Salz, Kupfer, Gips, Alabaster. - Tagungspublikation zum 17.  Treffen des Arbeitskreises Bergbaufolgeland- schaften vom 08.-09.04.2005 im Geopark Barbaros- sahöhle. Hannover, 225, pp. 12-17. Kupetz, M., 2008: Neue Vorstellungen zur Genese der Höhlen von Typ der Mansfeldischen Kalkschlot- ten.- In: Brust, M. K., Kupetz, M. & S. Schmiedel (eds.) Gips- und Anhydritkarst in der Mansfelder Mulde – Die Wimmelburger Schlotten, 23. Treffen Arbeitskreis Bergbaufolgen, Deutsche Gesellschaft für Geowissenschaften. 235, pp. 19-29. Kupetz, M., 2014.: Sulfatkarst am Kyffhäuser und die Genese der Barbarossahöhle.- In: Knolle, F. (eds.) Thüringen. Karst und Höhle 2011-2014, Verband der deutschen Höhlen- und Karstforscher, pp. 158- 163. Kupetz, M. & M. Brust, 2008: Historisches zum Begriff der Mansfeldischen Kalkschlotten sowie ein Beitrag zur nomenklatorischen Bestimmung dieses Höhlen- yps.- In: Brust, M. K., Kupetz, M. & S. Schmiedel (eds.) Gips- und Anhydritkarst in der Mansfelder Mulde – Die Wimmelburger Schlotten, 23. Treffen Arbeitskreis Bergbaufolgen, Deutsche Gesellschaft für Geowissenschaften. 235, pp. 61-74. Lebedev, A. & A. Lekhov, 1990: Dissolution kinetics of natural gypsum in water at 5-25 o C.- Geochem. Int., 27, 85-94. Lee, C.-H. & I. Farmer, 1993: Fluid flow in discontinuous rocks. Chapman & Hall, London. Li, J. & Z. Duan, 2011: A thermodynamic model for the prediction of phase equilibria and speciation in the H 2 O-CO 2 -NaCl-CaCO 3 -CaSO 4 system from 0 to 250 o C, 1 to 1000 bar with NaCl concentrations up to halite saturation.- Geochem. Cosmochem. Acta, 75, 4351-4376, https://doi.org/10.1016/j. gca.2011.05.019. Mehrbach, C., Culberson, C. H., Hawley, J. E. & R.M. Pyt- kowicx, 1973: Measurement of the apparent dis- sociation constants of carbonic acid in seawater at atmospheric pressure.- Limnology and Ocean- ography, 18(6), 897-907, https://doi.org/10.4319/ lo.1973.18.6.0897. ACTA CARSOLOGICA 48/2 – 2019 196 MODELLING SPELEOGENESIS IN SOLUBLE ROCKS: A CASE STUDY FROM THE PERMIAN ZECHSTEIN SEQUENCES EXPOSED ALONG THE SOUTHERN HARZ MOUNTAINS AND THE KYFFHÄUSER HILLS, GERMANY Miensopust, M., Hupfer, S., Kobe, M., Schneider-Löbens, C., Wadas, S. & C. Krawczyk, 2015: Geophysical in- vestigation of subrosion processes - an integrated approach.- Geophysical Research Abstracts, 18, pp. EGU2016-6799. Millero, F. J., 1979: The thermodynamics of the carbon- ate system in seawater.- Geochem. Cosmochem. Acta, 43, 1651-1661, https://doi.org/10.1016/0016- 7037(79)90184-4. Mucci, A., 1983: The solubility of calcite and aragonite in seawater at various salinities, temperatures, and one atmosphere total pressure.- Am. J. of Science, 283, 781-799, https://doi.org/10.2475/ajs.283.7.780 Parkhurst, D. & C. Appelo, 2013: Description of input and examples for PHREEQC version 3 - A computer pro- gram for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations.- U.S. Geological Survey Techniques and Methods 6, A43, pp. 497 [Online] Available from: http://pubs.usgs. gov/tm/06/a43. Plummer, L. N., Wigley, T. M. L. & D. L. Parkhurst, 1978: The kinetics of calcite dissolution in CO 2 -water systems at 5 o C to 60 o C and 0.0 to 1.0 atm CO 2 .- Am. J. Sci., 278, 179-216, https://doi.org/10.2475/ ajs.278.2.179. Scheffler, T. & T. Martienßen, 2013: Überwachungsmes- sungen zur Bestimmung der Deformationen von Kirchtürmen.- Schriftenreihe des Instituts für Mark- scheidewesen und Geodäsie der TU Bergakademie Freiberg, 1, 83-96. Schönenberg, R. & J. Neugebauer, 1987: Einführung in die Geologie Europas.- Rombach wissenschaft, 5 th edition. Freiburg, pp. 294. Schriel, W. & K. Bülow, 1926a: Geologische Karte von Preußen und benachbarten deutschen Ländern, Karte Bad Frankenhausen 4632.- Preußische geolo- gische Landesanstalt, Berlin, Germany. Schriel, W. & K. Bülow, 1926b: Geologische Karte von Preußen und benachbarten deutschen Ländern, Karte Kelbra 4532.- Preußische geologische Lande- sanstalt, Berlin, Germany. Serafeimidis, K. & G. Anagnostou, 2012: Simultaneous anhydrite dissolution and gypsum precipitation in a closed swelling rock system. In: Bobet, A. (eds.) 46th US Rock Mechanics / Geomechanics Sympo- sium, Chicago, pp. 3500. Spilker, M., 2016: Der Zechsteinkarst und sein Einfluss auf den Mansfelder Kupferschieferbergbau.- In: Grubenarchäologische Gesellschaft e.V. (eds.) 19. Internationaler Bergbau & Montanhistorik Work- shop Mansfeld-Südharz 2016. Clausthal-Zellerfeld, pp. 57-66. Svensson, U. & W. Dreybrodt, 1992: Dissolution kinet- ics of natural calcite minerals in CO 2 -water sys- tems approaching calcite equilibrium.- Chem. Geol., 100, 129-145, https://doi.org/10.1016/0009- 2541(92)90106-F. Völker, C. & R.  Völker, 1982: Die Elisabethschächter Schlotte.- Mitteilungen des Karstmuseums Uftrun- gen, 2. Völker, C. & R. Völker, 1991: Die Numburghöhle.- Mit- teilungen des Karstmuseums Heimkehle, 21. Völker, C. & R. Völker, 2014: Nach der grossen Entdeck- ung: Das Verschwinden der Numburghöhle im Karstwasser und das Vergessen der Ereignisse im Strudel von 1989.- In: Knolle, F. (eds.) Thüringen. Karst und Höhle 2011-2014, V erband der deutschen Höhlen- und Karstforscher, pp. 144-157. Wadas, S., Hupfer, S., Kobe, M., Miensopust, M. & C. Schneider-Löbens, 2017: Georisiko Erdfall: Und plötzlich war da ein Loch - Die Nachwuchsgruppe Subrosion des LIAG stellt sich vor.- DGG Mitteilun- gen, 01/2017, 17-23. Wadas, S., Polom, U. & C. Krawczyk, 2016: High-reso- lution shear wave reflection seismics as tool to im- age near-surface subrosion structures - a case study in Bad Frankenhausen, Germany.- Solid Earth, 7, 1491-1508, https://doi.org/10.5194/se-7-1491-2016. Weiss, R. F., 1974: Carbon dioxide in water and seawater: the solubility of a non-ideal gas.- Marine Chem- istry, 2, 203-215, https://doi.org/10.1016/0304- 4203(74)90015-2. Wessel, P . & W . H. F . Smith, 1998: New, improved version of generic mapping tools released.- EOS, 79, 579, https://doi.org/ 10.1029/98EO00426. Wessel, P., Smith, W. H. F., Scharroo, R., Luis, J.  F. & F. Wobbe, 2013: Generic Mapping Tools: Improved version released.- EOS, 94, 409-410, https://doi. org/10.1002/2013EO450001. Wissbrun, K., French, D. & P . Patterson Jr., 1954. The true ionization constant of carbonic acid in aqueous so- lution from 5 to 45 o .- J. Phys. Chem., 58(9), 693-695, https://doi.org/10.1021/j150519a004. Wunderlich, J., 2014: Zur Geologie des Kyffhäusers .- In: Knolle, F. (eds.) Thüringen. Karst und Höhle 2011- 2014, Verband der deutschen Höhlen- und Karst- forscher, pp. 138-143. Xu, P. & B. Yu, 2008: Developing a new form of perme- ability and Kozeny-Carman constant for homo- geneous porous media by means of fractal geom- etry.- Adv. Water Res., 31, 74-81, Developing a new form of permeability and Kozeny-Carman con- stant for homogeneous porous media by means of fractal geometry, https://doi.org/10.1016/j.advwa- tres.2007.06.003. ACTA CARSOLOGICA 48/2 – 2019 197