Salt dissolution potential estimated from two-dimensional vertical thermohaline flow and transport modeling along a Transylvanian salt diapir, Romania

The presented study explores the salt dissolution potential of fluids around a salt diapir in the Transylvanian town of Ocna Mures, Romania, which is facing land-collapse hazards related to salt mining activities. A structural three-dimensional (3D) model of the salt diapir, the adjacent basin sediments, and the mining galleries was developed based on existing maps, borehole data, own field observations, and geological publications of the Transylvanian Basin. The salt dissolution potential was estimated from 2D vertical thermohaline flow and transport model scenarios along the southeastern flank of the diapir. Results showed that the following factors increase the salt dissolution capacity along the upper 180 m of the diapir: (1) the presence of more permeable Quaternary alluvial sediments in connection with a fault zone of higher permeability along the diapir, and (2) the presence of more permeable sandstone units within the Miocene sediments in the east of the diapir, which provide freshwater access to the upper parts of the diapir. Thermohaline simulation with viscosity variation of the fluid, instead of a constant viscosity, influences the resulting salt fluxes by up to 50% within studied temperature ranges of 10–60 °C in the model domain. The range of theoretical dissolution rates along the upper 180 m of the diapir supports the hypothesis that cavern collapse is more likely to occur where cavern side walls have already been mined to almost no remaining side walls of rock salt, which is the case in the southeastern part of the diapir.


Introduction
Salt diapirs and the surrounding sediments are often involved in a variety of human activities such as salt mining, exploration and storage of hydrocarbons (e.g. Beckman and Williamson 1990;Hudec and Jackson 2007), and also repositories for radioactive waste material (Ludwig et al. 2001). The presence of highly soluble evaporitic rocks, a complex tectonic setting related to salt diapirsm, and human activities can lead to significant environmental problems, e.g. land subsidence, sinkhole development, salt cavern collapse, and contamination of water resources with brines.
The dissolution process of evaporites is closely linked to (1) the geological environment of the evaporitic deposits, (2) both the ability of freshwater to access the evaporites and the potential of brine to escape the system along preferential pathways, and (3) the presence of variable-density flow to accelerate the process (Johnson 1981). The structure of preferential pathways is inherited from geological structures such as discontinuities in sediments next to the salt dome, tectonic lineaments, which may influence permeability within fault zones, and the development of karst features within the evaporites. Circulation of fluids around a salt diapir plays a crucial role in the dissolution process, and the flow pattern itself is strongly influenced by the density variation of the fluids due to gradients of salt concentrations and temperature, respectively. The buoyancy ratio, which represents the relative effects of concentration and temperature on groundwater density, strongly controls whether the groundwater flow occurs either up or down along the flanks of salt diapirs (Evans and Nunn 1989;Evans et al. 1991). Evans et al. (1991) studied the effects of thermal and concentration gradients on two hypothetical salt domes. They found that if ambient formation-water salinities are high (greater than 130 g/l, resulting in a lower saltwater gradient), thermal effects become more dominant, and fluids can rise up along salt domes. Dissolved salt concentrations near the flank of salt domes can reach up to 300 g/l total dissolved solids (Ranganathan and Hanor 1989), which corresponds to a brine density of about 1,200 kg/m 3 (e.g. Van Duijn et al. 1988). Due to large variations in salinities and temperature, Magri et al. (2009) suggest to consider variable fluid viscosity in modeling studies around a salt dome in Schleswig-Holstein, Germany. They showed that using a constant fluid viscosity in the studied large-scale setting leads to overestimated salinity gradients.
A comparable small number of studies have demonstrated the importance of hydrogeological heterogeneity for variable-density flow around halites. Ludwig et al. (2001) studied the implications of large-scale heterogeneities reflecting differences in lithology at the potential site of radioactive waste at Gorleben (Germany) with two-dimensional/three-dimensional (2D/3D) groundwater models. They found a sensitivity of the saltwater/ freshwater system around the Gorleben salt dome to hydrogeological structures, which display a typical subvertical layering near the dome. Magri et al. (2009) studied the impact of more permeable transition zones in between a salt dome and adjacent units with coupled fluid-transport models. They found that vertical transition zones significantly influence the flow field, the distribution of heat, and the formation of highly concentrated saline plumes. In more recent studies on the effects of subsurface parameters on dissolution of halites, Zechner et al. (2011) andZidane et al. (2014a, b) showed that the presence of conductive vertical zones, which include a typical uncertainty related to the hydraulic characteristics of fault zones, plays an important role on variable-density flow of groundwater, transport of solutes, and dissolution of layered halite. Jamshidzadeh et al. (2015) studied density-driven flow around the Napoleonville Salt Dome in Louisiana, USA, and investigated the relationship of changes of groundwater flow, salinity, and heat transport due to a mining cavern failure at depths of more than 1 km. The cavern failure created a sinkhole of about 100 m diameter and more than 100 m depth overnight in August 2012. Jamshidzadeh et al. (2015) also investigated the effect of Boverpressured^fluids (i.e. pore fluid pressure is higher than corresponding hydrostatic pressure) on flow patterns in heterogeneous sediments around the diapir. Causes for overpressured fluids and their effects on flow patterns have been previously investigated, e.g. by Hanor and Bailey (1983), Nunn (1996) and McPherson and Garven (1999).
The presented study explores the salt dissolution potential of fluids around a salt dome in the Transylvanian town of Ocna Mures, which is facing land collapse hazards due to both the vulnerability of the near-surface rock salt, and its mining history. Rock salt in Ocna Mures has been explored since the Roman ages in open excavations, and up to the twentieth century in galleries and with solution mining. Most recently, in 2010 a sudden collapse in the adjacent Quaternary unconsolidated sediments led to the formation of a 70-90-m-wide salt lake with a depth of up to 23 m.
A joint Romanian-Swiss research project has developed geological and hydrogeological models in order to improve knowledge on possible hazards related to uncontrolled salt dissolution. One key aspect which is presented in this paper is the mechanism of density-driven thermohaline flow along the salt dome. Since the potential of salt dissolution is controlled by the concentration gradients and fluxes near vulnerable areas of the salt dome, specific attention has been given to hydrogeological boundary conditions: the hydraulic role of possible higher permeable structural fault zones along the salt dome, and the potential access to the salt dome of overpressured subsaturated groundwater within regional scale sandstone layers. A structural 3D model of the salt diapir, the adjacent basin sediments, and the mining galleries was developed based on existing maps, borehole data, own field observations, and geological publications of the Transylvanian Basin. The 3D model formed the base for 2D vertical density-dependent models of thermohaline flow and transport along the diapir. In combination with the 3D models of the mining structures within the salt dome, a first identification of possible areas with increased hazard is possible.

Three-dimensional model of hydrogeological and mining structures
The town of Ocna Mures, which is built around the salt diapir, is situated near to the western border of the Transylvanian Basin (Fig. 1). The basin consists of mainly Late Cretaceous to Late Miocene sediments (e.g. Krézsek et al. 2010;Tiliţă et al. 2013;Fig. 2). The salt diapir at Ocna Mures penetrated its cover of mainly shallow marine sediments and distinct tuff layers and reached a near surface position. The Quaternary cover is dominated by alluvial sediments of the post glacial river systems including the River Mures. The diapir itself evolved from a salt layer of early Badenian age (e.g. Ciulavu et al. 2000). Boreholes documented a base of the salt diapir in a depth between 1,000 and 1,700 m. The sedimentary rocks adjacent to the salt diapir are of middle to late Miocene age. They consist mainly of marls, sandstones, claystones and tuff units. While at Ocna Mures marls are dominating they are successively influenced towards the east by an increasing sand content forming intercalations of sandstone in the Middle Miocene (Fig. 2). The lithostratigraphic column of the Transylvanian Basin describes four sandstone units in between the distinct tuff layers BHadareni^on top, and BBorsa-Apahida, Iclod, Turda^at the base (Tiliţă et al. 2013). Freshwater springs have been described within sandstone layers east of Ocna Mures, indicating the possible presence of smaller sandstone aquifers within the Miocene marls (Posepny 1867).
The Ocna Mureş salt diapir is, unlike many other salt domes (e.g. Murray 1966), also characterized by an almost complete lack of a caprock cover. This could possibly be explained by the very high content of NaCl within the salt dome (over 98%; Stoica and Gherasie 1981), because caprocks are often an accumulation of less soluble minerals contained in the salt. The lack of less soluble minerals such as calcite, anhydrite/gypsum and clay would also result in a higher vulnerability of the diapir flanks to salt dissolution from subsaturated fluids. Huismans et al. (1997) acknowledged a compressional tectonic origin of salt diapirs structures in the Transylvanian Basin. However, no structures have been found in the field to document the stress field responsible for the diapirs. Distinct faults are not known in the area of Ocna Mures, but interpretation of seismic sections (Krezsek and Bally 2006;Onac 2010) crossing the Ocna Mures-Turda anticline indicate an eastward dipping thrust at which the diapir was possibly initiated, and which could border the Ocna Mures diapir along its eastern flank. A possible thrust of Miocene age, combined with a diapir rise may have led to a geomechanical weakening of sediments adjacent to the salt diapir. A possible increase of effective vertical permeability due to vertical faulting near rising salt domes together with groundwater being channelled along such faults has been discussed also by Evans and Nunn (1989). Therefore, the conceptual hydrogeological model used for the thermohaline simulations assumes for most scenarios a fault zone, which is hydraulically more permeable along the eastern border of the diapir.
The developed geological 3D model of the salt dome area includes all relevant available geological information in the area such as existing maps and cross-sections (Posepny 1867), borehole data, and geological publications of thẽ Based on topographical maps and on 24 2D cross sections, the six mine structures Stefania, located in the southeastern of the Ocna Mures diapir, Ferdinand, Francis and Josif (southwestern part of the diapir), 6th Martie (northwestern part of the diapir), and the Mine 1st Mai (northern part of the diapir) were constructed (Figs. 3 and 4). The mining chambers are mostly located in a zone in between 50 and 150 m depth below the top of the diapir. A horizontal model cross-section at 100 m depth shows that mine chambers Stefania in the southeastern part of the diapir were excavated the closest to the side walls of the diapir (Fig. 3b). Some of the southeastern edges of the chambers were even excavated into the marls outside of the diapir, and almost below the 2010 collapsed structure of Lake Plus. A NNW-ESE striking vertical cross-section through the diapir and the 2010 formed Lake Plus was extracted from the 3D model and used as a conceptual model for the development of the 2D thermohaline groundwater model scenarios (Fig. 4).

Two-dimensional vertical density dependent thermohaline flow and transport model east-southeast of the salt diapir
Aims of the 2D vertical density dependent flow and transport model were to investigate the possibility of density-driven thermohaline flow along subvertical fault zones of higher permeability next to the salt diapir, and the potential for subsaturated groundwater to dissolve the upper sides of the diapir. Due to the fact that former mining chambers in the SE (southeast) of the diapir were excavated very close to the boundary of the diapir, the potential for collapse of the side walls is higher in the SE of the diapir up to depths of more than 100 m. Different scenarios of 2D vertical flow and transport models were developed, which differ in (1) hydrogeological properties of permeable sandstone units in the Miocene sediments, of (2) the fault zone next to the diapir, and (3)  hydrostatic pressure imposed on permeable sandstone units and the fault zone. The scenarios were set up based on a WNW-ESE-striking 2D cross section through the diapir, the Lake Plus, and approximately 3,000 m into the Banta Hills, with a maximum depth of 1,300 m below sea level, and a maximum height of 425 m asl (Babove sea level^; see also Figs. 3b and 4). The spatial location of the 2D cross section was based on the following criteria: (1) to explore potential salt dissolution at the likely most vulnerable section of the salt dome flanks, which, according to the 2010 collapse leading to the formation of Lake Plus, is located in southeastern part of the dome and includes Lake Plus; (2) to include the modeled structures of the sandstone layers east of the dome due to their potential capacity to bring more Bcorrosive^freshwater close to the southeastern flanks of the salt dome; and (3) to extend the 2D section towards the ESE for up to 3 km in order to include the groundwater recharge area of the Banta Hills located SE of the dome. Groundwater recharge from comparably higher topographic elevations might explain overpressure, which creates increased flux of freshwater towards the flanks of the diapir. The limitation of 2D for the presented groundwater flow and transport scenarios is firstly due to the necessity to limit the number of unknown parameters in an area, where data on hydrogeological processes in the subsurface are very sparse. Secondly, simulation of density-driven thermohaline flow and transport on a larger regional scale with high density contrasts presents a challenging task even at the 2D scale. Finally, the chosen direct path from the flanks of the salt dome to the Banta Hills increases the possibility that the main direction of regional groundwater drainage towards the salt dome is parallel to the cross section.
The governing equations of variable-density fluid flow, mass, and heat transport were solved with the finite element software FEFLOW v6.2 (Diersch 2014). The 2D vertical models had between 20,000 and 25,000 triangular finite elements, depending on the scenario and the required mesh density for numerical stability: resulting sizes of triangle sides varied between less than 1 m to close to 50 m. Smaller elements with sides typically below 5 m were in and close to zones of higher hydraulic conductivities and, therefore, also were the highest density contrasts occurred (fault zone, sandstone layers, alluvials; Fig. 5a). The authors found that full upwinding was the only appropriate numerical scheme to face numerical oscillations and to ensure convergence of the solution. The typical drawback of the scheme, which is numerical dispersion (Diersch 2014), has been regarded as being of less importance, due to the fact that numerical dispersion is expected to be more important regarding the larger-sized elements, which are not located in areas where salt dissolution is simulated (boundary condition along upper parts of salt dome), and where the salt flux is higher (zones of higher hydraulic conductivities, see the preceding text). The density variation of the groundwater was coupled to the fluid concentration of rock salt NaCl, which resulted in a density range from 1,000 kg/m 3 (0 g/l NaCl) to 1,200 kg/m 3 (300 g/l NaCl). Viscosity variation of the fluid is calculated depending on simulated pressure, temperature and NaCl concentration with a relationship suggested by Diersch (2014), which is based on previous work from Lever and Jackson (1985), Hassanizadeh (1988), and Mercer and Pinder (1974).
The presented simulated scenarios were selected based on the idea to study the effect of unknown parameters such as hydraulic boundary conditions and hydraulic conductivities on salt dissolution at the vulnerable parts of the salt dome. Scenario 1 represents the model with the initial estimates of distribution of physical parameters and boundary conditions. Scenarios 2 and 3 differ from the initial scenario 1 with the variation of 1st type boundary condition for hydraulic head in order to explore the effect of under-, or overpressure within the fault zone and the sandstone layers. The next scenarios 4, 5 and 6 all focus on the effect of different (unknown) physical properties by lowering the hydraulic conductivities of the sandstone layers and the fault zone by up to two magnitudes. In addition, a scenario 1.1 with the same boundary conditions and model parameters as the initial scenario 1 was simulated with constant viscosity in order to estimate its effect on the simulated salt fluxes along the relevant upper 180 m of the diapir (see also section 'Results on near-surface scale'). Similarly, a scenario 1.2 with different geothermal parameters compared to the standard scenario 1 was simulated to study the estimated effect of (unknown) geothermal parameters on the resulting salt fluxes.
Model boundary conditions for fluid flow are prescribed head along the top of the model and follow the topographic elevation of the Banta Hills from 3 m below (topographic elevation) in valleys to 6 m below on hilltops, respectively (Figs. 4 and 5a). This approach was considered most appropriate considering the fact that no field data for groundwater recharge were available, and it follows basic approaches for regional groundwater flow (e.g. Hubbert 1940). Possible simulated flux errors due to the unknown effective piezometric head were considered of less relevance due the low permeability of the underlying marls (i.e. K = 1.0 × 10 −9 m/s), which results in small simulated fluxes. The four permeable sandstone units and the fault zone along the diapir were assumed to be 10 m wide. Both the four more permeable sandstone units and the fault zone have prescribed hydraulic heads at its east-southeastern end which are either (scenario 2) equal to the topographic elevation at that location, or are overpressured with prescribed heads of 350 m asl (scenario 1) and 500 masl (scenario 3). The 500-m asl overpressure from scenario 3 assumes a lateral hydraulic connection of the more permeable sandstone and fault zone units with the highest elevations of the Banta Hills, which in The highest elevation of the Banta Hills within the studies cross section reaches 425 m asl (Fig. 4). Boundary conditions for mass transport are a prescribed concentration of 300 g/l for the top of the halite of the diapir all along the bottom of the model, which corresponds to a maximum density of 1,200 kg/m 3 . Similarly, boundary conditions for heat transport are a prescribed temperature varying from 15 to 60°C along the diapir, which corresponds to an average geothermal gradient of 30°C/km depth increase.
Physical rock parameters were based on estimates for similar-type rocks, but could not be verified with field measurements at the site, except for conducted pump tests to estimate hydraulic conductivities of Quaternary sediments. For scenario 1, model parameters of hydraulic conductivities are K = 1.0 × 10 −5 m/s for the Quaternary alluvial foothill sediments (in the upper left model corner), K = 1.0 × 10 −7 m/s for the Miocene sandstone units, and the fault zone along the diapir, and K = 1.0 × 10 −9 m/s for the Miocene marls, respectively (Table 1). Uniform values have been assigned to the porosity (0.1), longitudinal dispersivity (10 m), transverse dispersivity (1 m), rock volumetric heat capacity (3.0 × 10 6 Jm −3 K −1 ), and rock thermal conductivity (2.0 Wm −1 K −1 ). Both rock volumetric heat capacity and thermal conductivity were lowered by 33% in scenario 1.2 in order to estimate these parameter sensitivities on simulated salt fluxes. Initial conditions for hydraulic heads, salt concentrations, and temperatures within the model domain were all set to 0. Simulation time was 10,000 years for all scenarios in order to reach close to steady-state conditions for thermohaline transport. To use the same geological 2D cross-section for longer simulations times has not been considered, also because geological boundary conditions around a potentially still-rising salt dome are nonstationary over longer time spans.

Results on model scale
Results on the model scale are presented for scenarios 1, 2, and 5, because these scenarios show exemplary differences in hydraulic heads and salt concentration (Fig. 5). Distributions of temperature do not vary significantly between scenarios and follow the imposed temperature gradient from 15°C at the top towards 60°C at the bottom of the model section, e.g. for scenario 1 (Fig. 5e). Simulated hydraulic heads after 10,000 years reflect the imposed boundary conditions along the topography, which reaches up to more than 400 m asl in the Banta Hills ( Fig. 5b-d). The increasing hydraulic head towards the model bottom in scenario 1, however, is due to the increasing fluid density, because the displayed equivalent freshwater heads as used by FEFLOW increase with the density ratio towards the bottom ( Fig. 5b; see also Diersch 2014). The reason for the increasing density ratio is visible by comparing it with the significantly increasing simulated salt concentration at the same location towards the bottom (Fig. 5f). Imposing lower hydraulic heads equivalent to the topographic elevation like in scenario 2 creates an underpressure in the   Tables 1  and 2 for definition of scenarios, and Table 2 for salt fluxes). Close-up of simulation results are shown in Fig. 6. Simulated hydraulic head values after 10,000 years are for scenario 1 (b), scenario 2 (c), and scenario 5 (d).
Simulated temperature values after 10,000 years are for scenario 1 (e); simulated salt concentration values after 10,000 years are for scenario 1 (f), scenario 2 (g), and scenario 5 (h) sandstone layers, and, moreover, in the lower part of the faultzone near the imposed boundary condition. Simulated concentrations after 10,000 years differ in between scenarios. As expected, highest concentrations up to saturation are simulated along the boundary condition along the flank of the diapir (Fig. 5f-h). Salinization within the model domain is dependent on the hydraulic boundary conditions imposed on the sandstone units of higher permeability: the lower imposed hydraulic head of scenario 2 leads to a Bcontamination^with denser saline water which propagates along the dip of the sandstone units, seeps towards the bottom of the model, and results in higher salinities in the lower model domain. Higher imposed hydraulic heads in the east of the sandstone units pushes freshwater along these units towards the upper part of the diapir reducing the fluid salinity and, therefore, increasing dissolution in these vulnerable parts (scenario 1, Figs. 5 and 6, and, even more so, scenario 3, Fig. 6). A scenario without any permeable sandstone units, and fault zone, respectively, results in a near-zero salinity in most of the model domain, except a very small and sharp concentration front along the flanks of the diapir (Fig. 5h).

Results on near-surface scale
The upper 180 m along the eastern flanks of the diapir represents the area suspected to be prone to collapse hazards. To understand possible causes of collapse, flow patterns and distributions of salt concentration are presented with close-ups of this area ( Fig. 6; see also the upper left model corner in Fig.  5a). The distribution of flow vectors at the end of the simulation offers a direct interpretation of flow patterns. For visibility reasons schematic flow vectors are used to summarize the nodal flow vectors simulated by FEFLOW after 10,000 years (white arrows in Fig. 6). The six scenarios are compared with the simulated salt flux from the prescribed concentration boundary condition along the upper (0-65 m), middle (65-130 m), and lower (130-180 m) part of the upper 180 m of the diapir, respectively. This corresponds also to the depths of former salt excavation chambers within the diapir (Figs. 3 and 4). Salt fluxes in g/day were integrated over each of the three parts at the end of the simulation time of 10,000 years for each of the scenarios. Integrated salt fluxes over the three parts Upper, Middle and Lower were sensitive to changes in the scenarios ( Table 2). The salt fluxes were also converted to potential subrosion rates (mm/year) along each the Upper, Middle, and Lower parts of the diapir based on a density of 2,170 kg/m 3 for halite ( Table 2). As expected, the simulated salt fluxes are strongly controlled by the exchange rate of groundwater along the boundary of the diapir, i.e. the access of freshwater towards the boundary of the diapir, and the transport of saturated saline water away from the boundary. Groundwater exchange is enhanced by the following factors: (1) increased hydraulic conductivity of the fault zone along the diapir, (2) a nearby connected larger freshwater reservoir like the more permeable Quaternary aquifer on top of the fault zone, and/or (3) larger pressure gradients due to boundary conditions representing increased regional groundwater recharge in connected sandstone layers. Salinity increase in the fault zone along the diapir typically leads to a downward flux due to the higher densities in almost all the simulated scenarios (Fig. 6). This results in a recharge of freshwater from the Quaternary aquifer on top, and, therefore, to increased salt fluxes, or subrosion in the Upper part along the diapir. Except if downward flux velocities are slowed down due to an upwelling of about equally pressured groundwater from the lateral sandstone units (scenario 1), or generally lower flux velocities due a decreased Fig. 6 Close-up on the concentration values towards the upper 200 m of the diapir for scenarios 1-6, the predominant flow directions within the more permeable units (if represented), the projected location of well 5200 (Fig. 3b), and the integrated boundary salt fluxes in g/day along the Upper, Middle, and Lower parts after 10,000 years simulation time. Location of close-up is shown in Fig. 5a Table 2   permeability in the fault zone (scenario 5). Larger overpressure in the sandstone units (scenario 3), or larger hydraulic conductivities in the upper part of the fault zone (scenario 6) lead to partly upward flow of freshwater in the fault zone, and, therefore to an increase of salt dissolution in that upper part, compared to the middle, or lower parts, respectively. Highest salt fluxes are simulated in the upper part when the sandstone units are underpressured, inducing a downward flow in both the fault zone and sandstones, and a compensation of more freshwater recharge from the Quaternary (scenario 2). Similarly, the model without more permeable sandstone units (scenario 4) creates an undisturbed downward flux in the fault zone inducing increased access of freshwater from the Quaternary.
Results of two additional scenarios were used to test model sensitivities to salt fluxes for simulations with constant viscosity (scenario 1.1) and varied geothermal parameters (scenario 1.2) compared to the standard scenario 1 (Table 2). Lowering geothermal parameters of volumetric heat capacity and thermal conductivity affects simulated salt fluxes not more than 12%. Using constant instead of variable viscosity in thermohaline transport simulation, however, affects salt fluxes in the upper, middle, and lower part by up to more than 50%, respectively.
Corresponding potential subrosion rates differ largely in between 0.1 mm/year (scenario 5; all parts) and 85 mm/year (scenario 2; upper part). All model scenarios, except scenario 1, indicate that the upper part is the most vulnerable to subrosion, and the lower part the less vulnerable (except scenario 5).

Comparison of simulations with measurements
The Ocna Mures study site offers almost no relevant data of hydraulic heads, salt concentration, or temperature distribution within the 2D vertical model domain for model calibration, or even model verification. However, measurements have been collected within the Quaternary aquifer, which is included in the upper few meters in the WNW corner of the 2D model scenarios. Because predominant groundwater flow in the Quaternary alluvials is horizontal, most of these measurements are influenced by flow directions not parallel to the 2D vertical plane of the model scenarios. Measured piezometric heads in well 5200, which is located less than 50 m southeast of saline Lake Stefania, and about 100 m south of saline Lake Plus (Fig.  3b), show values which are significantly higher than stage levels in both nearby lakes (Fig. 7), and also surrounding piezometers. Therefore, the authors assume that recharge to well 5200 is strongly influenced from saline vertical upwelling along the diapir, and within a possible fault zone. The projected model location of well 5200 is also approximately 300 m away from the nearest imposed 1st type boundary condition (of 257.00 m asl). Simulated piezometric heads and salt concentrations from model scenarios 1-6 and a model depth of 4 m are compared with measurements taken in 2014-2015 at a well depth of 4-5 m (Fig. 7). Simulated values of salt concentrations vary between 50 and 240 g/l (average measurements well 5200: 148 g/l), and simulated piezometric heads between 252.30 and almost 260 m asl (average measurements well 5200: 257.10 m asl). Scenario 6 with its more permeable upper fault zone (K = 1.10 −6 m/s) simulated almost the same values like the measurements. On the other hand, model scenarios with either strongly underpressured (scenario 2), or overpressured (scenario 3) sandstone layers under, or overestimate piezometric heads by up to 3 m, respectively. The model scenario with no permeable fault zone and sandstone units underestimates salt concentration by approx. 100 g/l, whereas model scenarios with overpressured sandstone layers (scenario 3), uniform fault zone hydraulic

Discussion
The discussion of the presented 2D thermohaline simulations focuses on the potential for subrosion of halite along the upper 180 m of the side walls, the limitations of the 2D simulations, and suggestions for further investigations. The upper 180 m coincide with the range of depths of the mining caverns within the diapir. 3D models of both the diapir structure and excavation cavern geometries have shown that the caverns in the SE part of the diapir were excavated up to the side walls of the diapir, and, at some locations, likely also into the Miocene sediments outside of the diapir (Figs. 3 and 4). Cohesion and stability of the Miocene sediments in the fault zone along the diapir could generally be reduced due to diapirsm and thrusting. The additional possibility of subrosion along the boundary of the diapir, and close to the remaining side walls of the caverns, might increase the hazard of chamber collapse. One geomechanical scenario for the Lake Plus collapse in 2010 could be an initial collapse of the nearest mining cavern within the diapir, the filling of the cavern with Miocene sediments from the fault zone, and the flushing out of the cavern brine into the newly forming salt lake at the surface. The presented 2D thermohaline scenarios give important indications of possible factors which contribute to hazards such as a cavern collapse. Thermohaline simulation with viscosity variation of the fluid appears to influence the resulting salt fluxes even within studied temperature ranges of 10-60°C in the model domain. The previously cited study from Magri et al. (2009), which called for including variable viscosity, has been studying temperature ranges in between 10 and 150°C. Simulated subrosion rates along the side walls of the diapir in the upper 130 m are typically several millimeters per year for most model scenarios. The already mentioned nearby mine structures Stefania in the southeastern part of the diapir are up to 150 years old. The most likely model scenario 6, according to the comparison with measurements in well 5200, reaches over the time span of 150 years a subrosion capacity of about 1-m side wall of the diapir for the upper 65 m, and 0.35 m in between depths of 65-130 m of the diapir side walls, respectively. Simulated subrosion rates of scenario 3, however, reach almost 13 m theoretical dissolution of diapir sidewalls in the upper 65 m in 150 years. This range of subrosion capacity confirms that cavern collapse due to lack of stability provided from rock salt side walls is more likely to occur where cavern side walls have been mined to almost no remaining side walls, which is the case in the southeastern part of the diapir near Lake Plus. The factors no, or small remaining cavern side walls, together with an active subrosion, increase the hazard for collapses along the studied southeastern flank of the diapir. The past land collapse from 2010, which formed the saline Lake Plus, has occurred in this area, which appears to be the most vulnerable to sudden land subsidence.
It is obvious that the limitation to 2D of the thermohaline simulations cannot account for hydraulic effects caused by the 3D structure of the diapir and a possible fault zone, heterogeneities of the more permeable sandstone layers within the Miocene sediments, and 3D regional groundwater flow conditions in the basin. Future developments of more complex geological, and hydrogeological models will rely on an additional exploration program of subsurface rock and fluid properties below the Quaternary sediments. Of primary focus would be the fault zone along the SE border of the diapir up to a depth of 180 m.

Conclusions
The presented study relies on results from 2D vertical thermohaline flow and transport simulations of saline groundwater along the east-southeastern flank of the Ocna Mures diapir. The simulations are based on a structural 3D model of the salt diapir, the adjacent basin sediments, and the mining galleries. Results of the simulations show that the following factors increase the salt dissolution capacity along the upper 180 m of the diapir: (1) the presence of more permeable Quaternary alluvial sediments in connection with a fault zone of higher permeability along the diapir, and (2) the presence of more permeable sandstone units within the Miocene sediments east of the diapir, which provide freshwater access to the upper parts of the diapir.
Salinity increase in the fault zone along the diapir typically leads to a downward flux due to the higher densities in almost all the simulated scenarios. This results in a recharge of freshwater from the Quaternary aquifer on top, and, therefore, to increased salt fluxes, or subrosion in the upper part along the diapir. Larger overpressure in the sandstone units, or larger hydraulic conductivities in the upper part of the fault zone lead to partly upward flow of freshwater in the fault zone, and, therefore also to a stronger increase of salt dissolution in this upper part.
The calculated theoretical dissolution rates along the upper 180 m of the diapir confirm that cavern collapse is more likely to occur where cavern side walls have been mined to almost no remaining side walls, which is the case in the southeastern part of the diapir. This might help to explain the sudden formation of a sinkhole filled with brine at the southeastern boundary of the diaper in 2010, and confirm the higher vulnerability to land collapses in this area.