Simulation of rock salt dissolution and its impact on land subsidence

Extensive land subsidence can occur due to subsurface dissolution of evaporites such as halite and gypsum. This paper explores techniques to simulate the salt dissolution forming an intrastratal karst, which is embedded in a sequence of carbonates, marls, anhydrite and gypsum. A numerical model is developed to simulate laminar flow in a subhorizontal void, which corresponds to an opening intrastratal karst. The numerical model is based on the laminar steady-state Stokes flow equation, and the advection dispersion transport equation coupled with the dissolution equation. The flow equation is solved using the nonconforming Crouzeix–Raviart (CR) finite element approximation for the Stokes equation. For the transport equation, a combination between discontinuous Galerkin method and multipoint flux approximation method is proposed. The numerical effect of the dissolution is considered by using a dynamic mesh variation that increases the size of the mesh based on the amount of dissolved salt. The numerical method is applied to a 2D geological cross section representing a Horst and Graben structure in the Tabular Jura of northwestern Switzerland. The model simulates salt dissolution within the geological section and predicts the amount of vertical dissolution as an indicator of potential subsidence that could occur. Simulation results showed that the highest dissolution amount is observed near the normal fault zones, and, therefore, the highest subsidence rates are expected above normal fault zones.


Introduction
Salt deposits (e.g.rock salt) are common in continental regions (Kozary et al., 1968).Different geo-mechanical problems, such as land subsidences and collapses, can be related to density-driven flow in evaporitic sedimentary rocks.Buried salt layers are dissolved and then removed due to the circulation of subsaturated groundwater (Quinlan et al., 1986;Frumkin, 2000;Hayek et al., 2012).Dissolution cavities that occur above or within the impermeable salt deposits lead to land subsidence or catastrophic collapses (Anderson and Kirkland, 1980;Reuter and Stoyan, 1993;Martinez et al., 1998;Frumkin and Raz, 2001).Moreover, salt dissolution also affects the water quality due to salinization and high mineralization (Johnson, 1981).The basic requirements for salt dissolution have been outlined and discussed by different authors, for example Johnson (1981), and Zechner et al. (2011).These requirements are (1) a deposit of salt or gypsum against which, or through which, water can flow, (2) a supply of water subsaturated with respect to NaCl or CaSO 4 , (3) an outlet where the resulting brine can escape, and (4) energy provided from hydrostatic head differences and/or density gradients, which causes groundwater flow through the system.Any combination of the four requirements can be induced and/or influenced by human activities such as mining, or groundwater pumping.Different aspects of salt dissolution can occur, as in some cases groundwater can interact with the salt rock from above (Johnson, 1981;Reuter and Stoyan, 1993), from the side (McManus and Hansor, 1993) or from below (Anderson and Kirkland, 1980).Gechter (2008) studied the geometry of dissolution cavities with rock salt dissolution experiments.His results showed that when freshwater gets in contact with an impermeable subhorizontal salt layer, a parallel void, or intrastratal karst formation develops in between the impermeable top layer and the salt layer.
Experiments have also shown that the spatial orientation of the geological boundaries between salt and non-soluble formations and the position of groundwater access play an important role for the shape of dissolution caverns (Gechter, 2008).Zechner et al. (2011) described the role of dip in geologic formations in a setting of horst and graben structures on salt dissolution: they found an increase of dissolution rates with increasing dips.Zidane et al. (2014) studied the influence of different subsurface parameters affecting the salt dissolution.The authors found that the number of faults with increased permeability within the normal fault zone is the most important factor that affects the dissolution compared to the other investigated parameters of aquifer thickness above the salt layer, a dynamic hydraulic conductivity of the aquifer above the salt layer, and varying boundary conditions due to pumping activities.They attributed the importance of fault heterogeneity on salt dissolution to the fact that the flow velocity increased at the top of the salt layer with the presence of up to six thinner faults compared to one large fault.
The hydraulic role of fracture systems, which represent pathways for freshwater access to salt formations embedded in less permeable formation, is still poorly understood.Therefore, the proposed models are extended to calculate intrastratal karst evolution as a function of hydraulically active fault systems.The calculations are based on simplified assumptions of fault zone geometries and intrastratal karst properties.The models consider an initial void, or fracture formation between an aquifer, where subsaturated groundwater can access the salt formation, and the salt formation itself.
One of the most commonly used techniques to model density-driven flow through fracture is based on the Navier-Stokes equation and coupled with the transport equation (Cardenas et al., 2007).When the flow is considered as laminar and steady, the inertial forces within the flow field are assumed to be negligible compared with the pressure and the viscous forces.In such cases, the flow equation is governed by the Stokes equation (Jäger and Mikelić, 2001).Due to the low velocity values recorded in previous studies within the aquifer above the salt layer (Zechner et al., 2011;Zidane et al., 2014), the flow is considered as laminar throughout the dissolution process and, hence, the Stokes equation is used to model the flow within the salt karst.The continuous development of the intrastratal karst creates an opening of the void at the top of the salt layer.To account for this phenomenon, a new formulation of the moving boundary condition is proposed to simulate the dissolution process with the laminar density-driven Stokes flow.The moving boundary condition is used to simulate the lowering of the void-salt interface due to dissolution at the interface.The moving boundary condition is based on a Dynamic Mesh Method (DMM) that adapts the size of the mesh with respect to the amount of displacement at the moving boundary.

Mathematical model
In the presented model only the liquid phase is simulated and the effect of dissolution is considered by increasing the size of the cells in the simulated domain with respect to the amount of dissolution that occurs.The idea behind the model is based on the fact that the dissolved volume of the solid gets continuously replaced by the fluid as the interface moves downward.Compared to different front capturing methods, the DMM does not require tracking the history of the front tracking points; only the coordinates of the boundary elements within the flow domain are modified which also reduces the high storage capacity required by most of the front tracking algorithms.The developed numerical model includes advanced approximations for both spatial and temporal discretization to reduce the computational needs and to maintain accuracy at the same time.The numerical code is used to simulate the density-driven flow problem including the dissolution equation.Flow and transport equations are coupled by the state equations linking density and viscosity variations to mass fraction.To achieve a high accuracy level for the spatial discretization of flow and transport, we used specifically suited methods.The flow through the opening void, or fracture is considered steady and laminar and the inertial forces in the flow field are assumed to be negligibly small compared with the viscous and pressure forces.Therefore, the free-flow is governed by the Stokes equation (Baumann and Oden, 1999;Saffman, 1971;Sanchez-Palencia, 1980;Kaviany, 1999;Jäger andMikelić, 2000, 2001).The Stokes flow is discretized using the Crouzeix-Raviart (CR) approximation, based on the nonconforming piecewise linear finite elements for the velocity and the piecewise constant finite elements for the pressure.This approximation provides locally mass conservative velocity which is an essential property for mass transport to avoid artificial mass sources and sinks.The transport equation is solved using a combination between the discontinuous Galerkin method (DG) and the multipoint flux approximation method (MPFA).The first one is used to discretize the advection part, while the second one is used for the discretization of the dispersion part.As the DG includes the advantages of both the finite element (FE) and the finite volume (FV) methods, it gives a robust and an accurate numerical scheme (Shuangzhang and Shahrouz, 2005).Considering the MPFA and as stated by Aavatsmark (2002), Aavatsmark et al. (1996), Edwards and Rogers (1998), Klausen and Russell (2004) and Wheeler and Yotov (2006), it is locally conservative and can handle irregular grids in anisotropic heterogeneous domains.On the other hand, both the MPFA and DG use the same type of unknowns that give the advantage of gathering it into one system matrix.So far, the combination DG-MPFA has shown to be a robust and accurate approach for modelling density-driven flow problems (Konz et al., 2008;Zechner et al., 2011;Zidane et al., 2012).Non-iterative time stepping is used in this work.The scheme is based on local truncation error control (Younes and Ackerer, 2010), which is able to increase the numerical accuracy and to reduce the computational cost at the same time.
Different models were developed to study the reactive transport problem.Graf and Therrien (2007) developed a 3-D numerical model to study heat and single-species reactive mass transport in a fractured porous media.The authors considered the fluid properties (density, viscosity) and the dissolution rate as a function of concentration and temperature.In our study we approximate the flow of a single species under 150 m depth with the specified boundary conditions (Fig. 1a) A. Zidane et al.: Simulation of rock salt dissolution and its impact on land subsidence as a function of concentration only.With the creeping flow within the intrastratal karst below the lower aquifer, we also consider the mass transfer coefficient as a constant during the simulation time.
The laminar groundwater flow and transport in the fracture can be described by the following system of equations (Flekkøy et al., 1996;Happel and Brenner, 1965;Landau and Lifshitz, 1987;Jäger and Mikelić, 2001;Cardenas et al., 2007).
The solute mass conservation equation is written in terms of mass fraction as follows (Boufadel, 2000;Cardenas et al., 2007;Younes et al., 2009;Zidane et al., 2012): State equations linking density and viscosity to mass fraction are The boundary conditions for the solute concentration at the moving boundary surface correspond to the first-order reaction equation where t is the time , ρ 1 and µ 1 the density and viscosity of the high density liquid (saltwater), respectively, ρ 0 and µ 0 the density and the viscosity of the liquid at zero concentration, and n is the unit outward normal to the boundary surface.

Spatial discretization of the flow equation
The system (Eqs.1-2) cannot be discretized with the same order for pressure and velocity approximations due to stability conditions.Otherwise some sort of stabilization needs to be added to the mixed formulation (Li and Chen, 2008).
To avoid these difficulties, we use the non-conforming Crouzeix-Raviart (CR) elements for the velocity approximation in combination with constant pressure per element, since they satisfy the Babuska-Brezzi condition (Brezzi and Fortin, 1991;Girault and Raviart, 1986;Gresho and Sani, 1998).This condition is central for ensuring that the final linear system to solve is non-singular (Langtangen et al., 2002).Moreover, the non-conforming Crouzeix-Raviart (CR) element has local mass conservation properties (Bruman and Hansbo, 2004) and leads to a relatively small number of unknowns due to the low-order shape functions.The CR element is used in many problems such as the Darcy-Stokes problem (Bruman and Hansbo, 2005), the Stokes problem (Crouzeix and Raviart, 1973) and the elasticity problem (Hansbo andLarson, 2002, 2003).The CR element gives a simple, stable and optimal order approximation of the Stokes equations (Arnold, 1993).In the following, we recall the main stages for the discretization of the Stokes equation with the CR triangular element.
With the non-conforming finite element method, the degrees of freedom for the velocity vector u are the two components (u i , v i ) of u at the mid-edge i. Inside the element E, we assume a linear variation of the velocity components (u E , v E ) where the interpolation function φ E i equals 1 on the mid-edge i and zero on the mid-edges j and k of E.
The variational formulation of the Stokes Eq. (1) using the test function φ i over the domain reads where ∇u is the gradient of the velocity vector u, and I is the 2 × 2 identity matrix.Using Green's formula gives The first integral contains boundary conditions.
Using Eq. ( 7) we obtain where x i =x j −x k and z i =z k −z j , z E and z i are respectively the z coordinate of the centre of E and of the midpoint of edge i, ρ E and p E are respectively the mean density and pressure over E.
The finite volume formulation of the continuity Eq. ( 2) over the element E reads: Using Eq. ( 6), it becomes The final system to solve for the flow is obtained by writing Eq. ( 9) for each edge (two equations per edge) and Eq. ( 13) for each element.

Spatial discretization of the transport equation
Standard numerical methods, such as FE or FV, usually generate solutions with numerical diffusion and/or non-physical oscillations when the advection part is dominant within the transport equation.The DG allows avoiding these oscillations (Siegel et al., 1997) since it provides a high-resolution scheme for advection.The local conservation of FV methods are maintained by the DG; in addition it allows higherorder approximations that could be used through a variational formulation in place of some hybridized difference or functional reconstruction (Kirby, 2000).The method was used on diffusion-advection problems in Baumann and Oden (1999), Cockburn et al. (1989), and Hugges et al. (2006), and multiple strategies have been used to adapt the DG method to elliptic problems (Aizinger et al., 2001;Arnold et al., 2002).More details of DG methods can be found in Arnold et al. (2002), and Cockburn et al. (1989) and Cockburn and Shu (1998).
Concerning the hyperbolic systems, the DG method has proven to be superior to the already existing FE methods (Arnold et al., 2002).The DG method is used to solve the advection equation and combined with the multipoint flux approximation (MPFA) method for the dispersion equation.The MPFA is locally conservative and handles general irregular grids on anisotropic heterogeneous domains.The MPFA method can be combined with the DG method without the time-splitting procedure (Younes and Ackerer, 2010).Because the MPFA and the DG use the same type of unknowns (average value per element), both discretizations can be gathered into one matrix system.
The spatial discretization of the DG-MPFA is given as follows: by substituting the mass conservation of the fluid in the transport equation, the transport Eq. ( 3) can then be written in the following mixed form: The dispersive flux u D is assumed to vary linearly inside the element E, therefore, where Q E D,i is the dispersive flux across the edge i, and Q s is the sink/source term.
We use the P 1 DG method where the approximate solution C h (x, t) is expressed with linear basis functions ϕ E i on each element E as follows: The three unknowns for each element are the average value of the mass fraction defined at the triangle centroid (x E , y E ) and its deviations in each space direction (Cockburn and Shu, 1998) with the corresponding interpolation functions: The variational formulation of Eq. ( 14) over the element E using ϕ E i as test functions gives where ∂E is the boundary edge of the element E, and C * the upstream mass fraction on ∂E.

Dissolution process
To model the salt dissolution process, we used a technique based on the variation of the size of the mesh.The simulated domain is only the water circulation area (i.e.flow channel), and as long as the dissolution occurs, the size of the mesh within this area will increase (Fig. 2).The variation of mesh size is directly related to the amount of dissolved salt.For an element E with area |E|, the dissolved mass within an interval of time dt is given as follows: where dm is the amount of the dissolved mass within an interval of time dt.Knowing the density definition of a certain amount of salt results in where ρ s is the salt density, dA•e is the volume of the dissolved salt, with dA the dissolved area and e the dissolved thickness.In the case of a 2-D dissolution process the volume is then reduced to the area dA.Using Eq. ( 19) in Eq. ( 20) The area of the dissolved salt can be approximated as follows (Fig. 2): where dh is the incremental height that is added to the edge at the salt boundary, and dy is the height of that edge.In the case of multiple edges at the salt boundary, and since the coordinate variation is related to the nodes, the amount of dissolved area at each edge is divided on the two corresponding nodes.In this case, the height dy used in Eq. ( 22) to deduce the amount of coordinate increase dh for each node is nothing but the sum of the halves of the two edges sharing the same node (Fig. 2).Consequently, the increase dh can be calculated for each node n at the salt boundary as follows: The increase dh is subsequently added to the horizontal coordinate of that node at each time step.

Coupling flow and transport equations
Numerical simulations of density-driven problems require excessive computational time due to the strong nonlinearities between the flow and the transport equations.In order to reduce the computational needs and maintain accuracy at the same time, a non-iterative time-stepping scheme based on local truncation error control is used in this study.Hirthe and Graf (2012) implemented the non-iterative second-order time-stepping scheme of Kavetski et al. (2002) on a Hydro-GeoSphere model to solve the low Rayleigh number-thermal Elder problem of free convection in porous media.In the presented work we use an adaptive time stepping based on the local truncation error of the concentration (Younes and Ackerer, 2010).The time-stepping procedure is shown as follows.
The local time truncation error of the concentration is evaluated using two approximations of adjacent order of accuracy.
The time step is accepted if the absolute error criterion is verified: If this criterion is met, the following time step is controlled by the temporal truncation error tolerance γ using If the error criterion is not satisfied, the current time step is repeated using the latest error estimate  often set equal to 2.0 and 0.1, respectively, and s = 0.9 is a safety factor (Sloan and Abbo, 1999).The temporal truncation error tolerance is set to γ = 0.01.

Field scale 2-D cross section
To study the influence of fault systems on intrastratal karst formation, a 2-D vertical model representing Mesozoic sediments structured in Horst and Graben structures typical for the Tabular Jura of northwestern Switzerland was set up (Zechner et al., 2011).The domain has the characteristics shown in (Table 1).The section is constituted by an upper aquifer and a lower aquifer that are connected via two fault zones at the WNW and ESE sides, respectively.The 2-D cross section with 1000 m length and 150 m depth is based on an existing three-dimensional geological model of 47 faults and 4 faulted horizons of the main aquifers and aquitards and a derived regional groundwater flow model (Spottke et al., 2005).A constant hydraulic head is imposed at the ESE side, and a pumping well with constant head is applied at the WNW side (Fig. 1a).The pumping rate is 0.15 m 3 s −1 , which is equivalent to a constant head of 251 m.In different related studies (Zechner et al., 2011;Zidane et al., 2014) the effect of the pumping rate and the well depth (into the aquifer) were evaluated in detail.The effect of the pumping in a 2-D domain for a well of length L (depth into the aquifer) is considered by dividing the pumping rate over the finite elements that intersect with the well within the upper aquifer.
The individual pumping rate of each finite element is then assigned to the sink/source term (Q s ) in the transport equation.A constant head is applied at the ESE boundary of the upper aquifer and a pumping well which creates a constant head is imposed at the WNW boundary.For the transport boundary conditions a zero constant relative concentration is imposed at the inlet ESE boundary, whereas the relative concentration is set to one at the top of the salt layer, that is the bottom of the simulated domain (Fig. 1a).The salt formation is located within the lowest layer of the lower aquifer, and therefore dissolution can only occur at the bottom of the section.To simplify the modelling approach, 1 mm of salt is considered as previously dissolved and the fracture evolution is considered within this initially developed void.The 1 mm thickness of dissolved salt is then removed from the final aperture width when calculating the vertical dissolution within the lower aquifer.The lowest layer of the salt formation is discretized by 937 triangular finite elements.The amount of dissolved salt and the corresponding displacement of the top of salt are calculated within each element of the salt formation over the 30 years of simulation.Three different scenarios were evaluated.In the first case, the fault thickness at the ESE and WNW sides is set to 10 m.In the second case, the fault thickness is increased to 40 m, and finally, in the third case, the 40 m fault is replaced with six thinner faults.The dissolution rate and corresponding vertical displacement of the top of salt are calculated within the three studied cases.

Test case 1
In this test case, the fault thickness at the ESE and WNW sides are set to 10 m.The boundary conditions are kept the same as discussed before.The dissolution rate is calculated over the lowest layer of the lower aquifer, and results show that the highest amount of dissolved salt is recorded near the fault zones (Fig. 3).With increasing distance from the fault zones the dissolved mass decreases sharply after less than 100 m distance from the fault zones (Fig. 3).The amount of vertical dissolution is calculated after 1 year and after 30 years of simulations (Fig. 4).Results show that the highest displacements are recorded near the fault zones where the subsidence reaches up to 30 mm (Fig. 4).The vertical dissolution reduces gradually with distance from the fault zones and remains at an almost constant and negligible rate between 0.6 and 0.7 mm within the lower aquifer (Fig. 4).The concentration distribution within the lower aquifer shows that near the fault zones the water is less saturated, which enhances the dissolution and therefore possible land subsidence (Fig. 5).With increasing distance from the fault zones water is completely saturated in the lower aquifer, which explains the fact that dissolution is negligible within the 100 to 900 m segment of the lower aquifer (Figs. 4 and 5).

Test case 2
For the second test case, the fault width at the ENE and WNW end of the section is increased to 40 m with the same boundary conditions as discussed before.With 40 m width, low concentration values are observed near the fault zones (Fig. 6), which indicates that high concentration gradients are expected near these areas.Compared to the previous test case, a different dissolution profile is observed within the lower aquifer in the case of 40 m wide fault zones.In this case the dissolution rate within the lower aquifer is no longer negligible.Indeed, a remarkable vertical dissolution is observed within the lower aquifer up to 400 m of distance from the WNW side within the lower aquifer (Fig. 7).Negligible vertical dissolution is recorded within the 400 to 600 m segment.The simulated dissolution increases to reach its maximum value at the ESE side of the domain (Fig. 7).The amount of vertical dissolution is similar to the dissolution profile.Important displacement values are recorded near the faults, and they remain at a significant rate within the lower aquifer to reach their minimum values within the 400 to 600 m segment of the section (Fig. 8).It is remarkable, however, that the maximum possible vertical dissolution with 40 m wide fault zones does not exceed the maximum values with 10 m fault thickness (0.03 m).Thus, with a wide fault the maximum possible vertical dissolution does not exceed the values of the previous test case, whereas the area that is affected by vertical dissolution might be much larger (Fig. 8).

Test case 3
In the last test case the 40 m wide fault was replaced with multiple thin faults at the ESE and WNW sides.The 40 m wide fault is replaced with six faults as follows:= -1 fault -1 × 40 m thickness -0 m gap -6 faults -6 × 2.5 m thickness -5 × 5 m gap.
Similar to previous test cases, the concentration values show that groundwater is fully saturated with NaCl far from the faults and low concentration values are recorded near the six thin faults at both the ESE and WNW sides (Fig. 9).A high dissolved mass is observed compared to previous test cases (Fig. 10): the dissolved mass reaches up to 200 kg as a maximum value and reduces gradually to reach negligible values within the 200 to 800 m segment within the lower aquifer.The possible vertical dissolution is 3 times higher compared to previous test cases (Fig. 11).In fact, a 0.1 m vertical displacement is observed at the ESE side of the aquifer, which gradually reduces with distance from the faults to reach negligible values within the 200 to 800 m segment, and then increases to 0.04-0.06m in the 0-100 m segment at the WNW side of the aquifer (Fig. 11).When comparing the dissolution profiles of the three different test cases it could be observed that the dissolution profile is symmetrical within the first test case, but not within the last two test cases (Figs. 5 and 9).This is due to the fact that with larger and/or multiple faults, the effect of the fixed hydraulic head boundary at the ESE side is more important than the effect of the pumping well at the WNW side in increasing the concentration gradient (i.e.reducing the NaCl saturation) above the salt layer in the lower aquifer.

Discussion
The presented approach simulates the Stokes flow in a subhorizontal intrastratal karst and uses the moving boundary method to account for salt dissolution-coupled void opening.In all the simulated test cases the flow directions show a similar pattern (Fig. 12).Due to the higher hydraulic head imposed at the ESE boundary, freshwater enters the upper aquifer from the ESE, and flows towards the well.Within the lower aquifer, however, flow direction is from WNW to ESE due to the gravitational forces of the denser, more saline groundwater.The mainly freshwater and partly saltwater mixture is extracted by the pumping well at the WNW side of the upper aquifer.The assumption that the flow remains laminar during all the simulated dissolution process appears to be justified due to the low velocity values also recorded in previous studies within the lower aquifer (Zechner et al., 2011;Zidane et al., 2014) where the highest values recorded were also in the range of 4-5 m d −1 .The void opening, which occurs mainly in the first 3 years of the simulation, leads to a further lowering of the flow velocities and therefore the flow remains laminar after the dissolution process within the lower aquifer; hence the Stokes equation is applied within the void throughout the entire simulation time.
In addition, a new formulation of the moving boundary condition is proposed to simulate the dissolution process with the laminar density-driven Stokes flow.The presented method allows for simulating void opening, which results in a vertical lowering of the top of the salt formation.This simulated lowering of the illustrated test cases does not automatically correspond to subsidence rates, which are observed at the ground surface.This is due to the geo-mechanical response of the overlying strata to the fracture opening (Fokker, 1995), which has not been investigated in this study.Nevertheless, the temporal and spatial distribution of the simulated displacement along the 2-D cross sections can be qualitatively compared to the shape and subsidence rate of the observed subsidence areas in northwestern Switzerland (Fig. 1a and b).From the six observed subsidence areas, three are likely not directly related to solution mining, but rather to planar subsurface dissolution of evaporites (west, north, central).Subsidence rates as estimated over an observation period from 1975 to 2000 are up to 0.012 m year −1 in the west and the north, respectively, and up to 0.12 m year −1 in the central area.These rates are comparable to the simulated vertical displacements over the first year: 0.01 to 0.02 m year −1 for the test cases with single fault zones, and up to 0.1 m year −1 for multiple faults connecting the upper aquifer with the fracture above the salt top.At the time t = 0, the aquifer is initially filled with freshwater, hence the concentration is set to zero in the entire domain, including the lower aquifer.Hence, the high concentration gradient (∇C) at the early time steps induces the high dissolution rate at the beginning of the simulation.We should point out that this study is a simplification of the real setup that could be  observed in field scale.The vertical dissolution rate in the simulated test cases is not constant over the 30 years of simulation.About two-thirds of dissolution occurs in the first year, and after 3 years the simulated dissolution does not increase anymore.The reason for the decrease of dissolution over time is that the system reaches steady-state flow conditions and, therefore, no subsaturated groundwater reaches the top of the salt anymore.The observed subsidence yearly rates, however, are averaged over a period of 25 years, and add up to absolute subsidences which are much higher than the amounts of simulated vertical dissolution.This suggests that the hydrogeological properties and boundary conditions at the studied site are likely more transient and, therefore, constantly open new pathways for subsaturated water to access the salt top, and saturated water to leave the salt top.
The observed diameters of the affected areas are more than 1000 m in the west, 300-500 m in the north, and less than 200m in the central area.Different lengths of dissolutionaffected parts can also be observed along the 2-D cross section: the wider the hydraulically active fault zone, the higher the influence of dissolution with increasing distance from the fault zone.The simulations point towards the possibility of the western area being characterized by wider fault zones, whereas the "small-area" and "high-rate" subsidence type in the central area might be related to a complex fault zone of several thinner fault zones.The mechanical behaviour of the salt overburden, such as interbedded marls and clay stones, carbonates and interbedded marls and evaporite-bearing sequences, could show a complex behaviour such as brittle or ductile in space and time (Geluk et al., 2007).There are also indications that the variation of the thickness in the rock salt layer in the Tabular Jura is not only of sedimentary origin, but likely to be the result of normal-and strike-slip faulting (Spottke et al., 2005), including extension and subsequent faulting which could lead to weakening in the sedimentary cover and, therefore, enhance karstification, particularly near fault zones.The effect of fault zones on subsidence was also documented by Kruse (1999): they observed an increase of subsidence in numerical experiments with fault zones.
Although the present work represents a simplification of a complex weathered fault system, the results constitute a base for rock mechanical stability calculations.In order to compare future rock mechanical stability calculations with subsidence rates, a continuous high-resolution geodetic monitoring in time and in space (i.e. also in between "well-known" subsidence areas) will allow further understanding of the mechanisms of subsurface intrastratal salt karst development in populated areas.The equation of dissolution is then added to the model and the dissolution process is simulated by using a dynamic mesh routine that adapts the size of the mesh with respect to the amount of dissolution that occurs.
A 2-D cross section with 1000 m length and 150 m depth based on field measurement, a conceptualization of the setting in the Tabular Jura near Basel, is used to simulate the salt dissolution within the underlying halite formation beneath the lower aquifer.The developed numerical model is used to calculate the amount of increase that should be added to the finite elements on top of the salt formation due to the dissolution process.
Results show that highest void openings are recorded near the fault zones and decrease gradually to a minimal value towards the central part of the lower aquifer.Three test cases were studied.In the first test case a single 10 m wide fault zone at both the ESE and WNW sides of the 2-D cross section was used.In the second test case the 10 m faults at the ESE and WNW sides were replaced with 40 m wide faults.In the last test case the 40 m wide fault is replaced with six thin faults.The study shows that with large faults (40 m) the total amount of dissolved salt is higher than with the 10 m wide faults.The only difference between the first and the second test cases regarding the expected dissolution rates is the relatively high displacement rates (0.005-0.01 m) along the salt top far from the fault zones.However, with multiple faults as simulated in the third test case, very high amounts of displacements (i.e.vertical dissolution) are observed at the ESE (0.1 m) and WNW (0.02-0.04 m) sides of the domain.The present study shows that the velocity values do not exceed the 5-6 m d −1 within the karst formation, which corresponds to low Reynolds number (Re < 100) in both cases, hence, the Stokes equation is considered applicable to model the laminar flow within the salt karst.
The paper explored numerical techniques to simulate the salt dissolution and represent an important boundary condition when calculating the linear (elastic) and non-linear mechanical behaviour of the overlying geological formations such as carbonates and salty-marls, respectively clay stones to the development of an intrastratal karst.

Figure 1 .
Figure 1.(a) Domain and boundary conditions of the 2-D cross section, and conceptual model for simulation of the dissolution process within the (initial) fracture.(b) Map of the study area in northwestern Switzerland with location of 2-D cross section.Interpolation of subsidence rates are between 1975 and 2000 (averaged over 25 years), with focus on central (lower left panel) and northern (lower right panel) subsidence areas, respectively.
et al.: Simulation of rock salt dissolution and its impact on land subsidence

Figure 2 .
Figure 2. Area and relative height dh of the dissolved salt at the boundary between the fracture and the salt layer.

Figure 3 .
Figure 3. Dissolved salt over each element after 3 years of simulation with 10 m fault thickness.

Figure 4 .
Figure 4. Amount of vertical dissolution over each element after 1 and 3 years of simulation with 10 m fault thickness.

Figure 5 .
Figure 5. NaCl concentration distribution over the domain after 3 years of simulation with 10 m fault thickness.

Figure 6 .
Figure 6.NaCl concentration distribution over the domain with 40 m wide fault zone after 3 years of simulations.

Figure 7 .
Figure 7. Dissolved salt over each element after 3 years of simulation with 40 m wide fault zone.

Figure 8 .
Figure 8. Amount of vertical dissolution over each element after 6 months and 3 years of simulation with 40 m wide fault zone.

Figure 10 .
Figure 10.Dissolved salt over each element after 3 years of simulation with 6 each of 15 m wide faults.

Figure 11 .
Figure 11.Amount of vertical dissolution over each element after 6 months and 3 years of simulation with 6 thin faults.

Figure 12 .
Figure 12.Flow vectors within the 2-D cross section for test case 1 after 30 years.

Zidane et al.: Simulation of rock salt dissolution and its impact on land subsidence 2181
It is eliminated in the case of a free-flow boundary or in the case of an Hydrol.Earth Syst.Sci., 18, 2177-2189, 2014 www.hydrol-earth-syst-sci.net/18/2177/2014/ A.interior edge i.In this last case, Eq. (8) becomes

Table 1 .
Values of physical and geometrical parameters of the problem.