A Meshfree Generalized Finite Difference Method for Solution Mining Processes

Experimental and field investigations for solution mining processes have improved intensely in recent years. Due to today's computing capacities, three-dimensional simulations of potential salt solution caverns can further enhance the understanding of these processes. They serve as a"virtual prototype"of a projected site and support planning in reasonable time. In this contribution, we present a meshfree Generalized Finite Difference Method (GFDM) based on a cloud of numerical points that is able to simulate solution mining processes on microscopic as well as macroscopic scales, which differ significantly in both the spatial and temporal scale. Focusing on anticipated industrial requirements, Lagrangian and Eulerian formulations including an Arbitrary Lagrangian-Eulerian (ALE) approach are considered.


INTRODUCTION
The basic motivation of this research is to provide a method that is able to simulate the long-term development of a salt cavern during a double-well solution mining process. Solution mining is used to extract underground water-soluble minerals such as salt and potash. A double-well convection process has been a preferred choice for solution mining due to it's large recovery rate [2,33]. As the name suggests, this involves the use of two boreholes or wells for the extraction process: an injection well, and a recovery or extraction well. For the extraction of salt, fresh water is pumped into a salt deposit through the first "injection" well. Salt present in the cavern dissolves in the water to produce a saturated brine solution. This is then extracted at the second "extraction" well. A schematic of this process is shown in Fig. 1. The main direction of dissolution is vertical, which is controlled by alternate lifting of the injection and the extraction well.
In this work, we focus on modeling of the fluid flow involved in such a double-well solution mining procedure, including the formation of the salt-water solution. An essential aspect of this is to accurately model the long-term geometrical evolution of the salt cavern. This is needed to steer the actual process of solution mining, in terms of, for example, determining when and at what rate the injection and extraction wells are raised. However, numerically modeling this is very challenging, as it is a highly dynamic three-dimensional process involving different spatial and temporal scales. Over the time scale of several years, as salt in the cavern dissolves in the water, the cavern starts to erode, causing significant deformations in it's overall shape. However, the dissolution process relies on a smaller time scale of several minutes. On the spatial scale, the former involves the modeling of the entire salt cavern, while the latter is more localized and is relevant near the cavern walls. In the present work, we model both these processes in separate simulation setups. A macroscopic simulation is carried out to model the evolution of the cavern over many years. This is done on actual salt cavern geometries. The computation Figure 1. Schematic of double-well solution mining (adapted from [25]). of the diffusion rate of the salt (and related minerals) to be used in these macroscopic simulations are done in a separate simulation, in the so-called microscopic setup. This involves simulations over the smaller time scale of a few minutes, and over representative geometries several orders of magnitude smaller than the size of the salt cavern in the macroscopic simulations.
Over the past few decades, meshfree simulation methods have emerged as an alternative to the conventionally used mesh-based simulation procedures, especially in the context of modeling fluid flow. The advantages of meshfree methods are most notable for applications with complex domains, or those with moving geometry parts, free surfaces, phase boundaries, or large deformations. While modeling each of the latter cases with mesh-based methods, the highly dynamic nature of the simulations often requires an expensive global remeshing procedure. On the other hand, moving Lagrangian and semi-Lagrangian procedures fit in naturally with meshfree methods, making the simulation of dynamic geometries or phase boundaries a lot easier. In the application at hand, the modeling of the changing domain during the long-term evolution of the salt cavern falls into this category. We thus choose a meshfree approach.
In this paper, we use a meshfree Generalized Finite Difference Method (GFDM) [5,7,12,20] based on a cloud of numerical points. This method has already been successfully applied in various CFD and continuum mechanics applications. Prominent examples include water crossing of cars, water turbines, hydraulic valves, soil mechanics, metal cutting, and mold filling [10,14,21,23,31,32]. The methods presented here are part of the in-house developed software MESHFREE † , which combines the advantages of GFDM and the fast linear solvers of SAMG [22].
We start by using a Lagrangian formulation where the discretizing point cloud moves according to the flow velocity [10,13]. This results in an accurate and natural transport of physical information. The basic physical model consists of the conservation equations for mass, momentum, and energy. For solution mining processes, we extend it by the standard k-ε turbulence model and equations for the concentration of the occurring species (see Sect. 2). The GFDM specific numerics are presented in Sect. 3 with special focus on the Lagrangian and Eulerian formulations. Microscopic simulations are presented in Sect. 4, and are used to determine the necessary effective model parameters of the macroscopic problem which follows. For macroscopic simulations, the Lagrangian formulation leads to a significant restriction of the time step size due to the explicit movement of the point cloud. To enable simulations in reasonable time, an Eulerian formulation should be preferred in this context. Here, the point cloud is fixed and convective terms represent the transport of physical information. The movement of the boundary of the salt cavern is based on the solution rate of the salt in the flowing water. To accurately handle this moving boundary, interior points close to the boundary are subject to an ALE-approach (Arbitrary Lagrangian-Eulerian) according to [8]. This procedure gives rise to covering the complete life cycle of a salt cavern -several decades -by a meshfree simulation. In Sect. 5, we demonstrate the advantages of the Eulerian formulation for a simplified macroscopic example of a double-well solution mining process, followed by concluding remarks in Sect. 6.

PHYSICAL MODEL
In this section, we describe the basic physical flow model and its extensions for modeling solution mining processes, in both the macroscopic as well as microscopic simulations. Specific models needed for the density, viscosity, and heat capacity of a solution are also discussed.

Basic equations
The basic underlying physical model is given by the conservation equations of mass, momentum, and energy in Lagrangian formulation.
In general, the stress tensor is split into its viscous and solid parts by S = S visc + S solid [10,13]. For the present application, the stress tensor is purely viscous, S solid = 0. The viscous part is defined by where I ∈ R 3×3 is the identity. To incorporate turbulent effects, the standard k-ε turbulence model [18] is considered for turbulent kinetic energy k and turbulent dissipation ε where η is the laminar viscosity, and η turb = ρ ·C η · k 2 ε is the turbulent viscosity. Fluctuating dilatation and source terms are omitted [18]. The turbulent production rate is defined by P pr = η turb · ∇v T 2 M with von Mises matrix norm · M . The turbulent buoyancy is given by Pr turb · ∂ ρ ∂ T ·(g·∇T ). For this model, well-established values for the constants are used σ k = 1.0, σ ε = 1.3, C 1ε = 1.44, C 2ε = 1.92, C 3ε = −0.33, C η = 0.09, and turbulent Prandtl number Pr turb = 0.85. Furthermore, a logarithmic wall function is used in the vicinity of walls.
In order to simulate solution mining processes, the basic model above is extended by convectiondiffusion equations to represent the different minerals or species present in the salt mixture. For the concentration c i of species i = 1, . . . , N with effective diffusion coefficient D i,eff , we have In the Eulerian formulation, the material derivative is replaced by its definition, i.e. d dt = ∂ ∂t + v T ∇.

Modeling density, viscosity, and heat capacity
Consider the general form of the equation of state where density depends on the temperature and each of the concentrations. Based on the formulation in [16,17], the density of a solution of N species in water is given by where w H 2 O and w i are the mass fraction of water and species i, respectively. Additionally, w H 2 O + ∑ N i=1 w i = 1 has to be satisfied. The density of water is determined by the non-linear relation with A 1 , . . . , A 7 and C * 0 , . . . ,C * 4 defined according to [17]. The apparent density of species i is given by The mass fractions w i are defined by the concentrations c i as The viscosity of the solution, η sol (T, c 1 , . . . , c N ), and its heat capacity c v,sol (T, c 1 , . . . , c N ) are modeled in a similar manner. For the viscosity of a solution of N species in water, we use a modified version of the Arrhenius equation where the viscosity of water depends on the temperature as Furthermore, the viscosity of species i is given by with constants V * 1 , . . . ,V * 6 according to [15].
Furthermore, the heat capacity of species i is modeled by where and constants B * 1 , . . . , B * 6 are according to [16]. We use quadratic interpolation (extrapolation) for the definition of the heat capacity of water. Assume given temperatures T 1 , T 2 , T 3 , with T 2 = T 1 + ∆T and T 3 = T 2 + ∆T (∆T > 0). The corresponding heat capacities of water c v,T 1 , c v,T 2 , and c v,T 3 are also assumed given. Then, the heat capacity of water at arbitrary temperature T is determined by T 1 , T 2 , T 3 are chosen adaptively with ∆T = 5 • C, depending on the value of T . The range of c v,T k values, k = 1, 2, 3, are taken from [16], which provides the values between 0 • C and 95 • C. We restrict the study in this paper to sodium chloride as the species of interest. All the model constants mentioned in this section for sodium chloride are summarized in Table I.

Point cloud preliminaries
In the GFDM approach, the computational domain is discretized by a cloud of numerical points. The point cloud is composed of NP = NP(t) number of points, which includes points in the interior of the domain, and those at the boundary. The initial seeding of these point clouds is done by a meshfree advancing front technique, details of which can be found in [19,24]. The density of the point cloud is given by a sufficiently smooth function h = h(x,t), the so-called interaction radius or smoothing length. Thus, h prescribes the resolution of the point cloud. It is also used to define the neighborhood of each point. For a point x j in the point cloud, all approximations are performed using only nearby points within a distance h from it. This set of nearby points is referred to as the neighborhood or support of x j , and is denoted by To ensure a sufficient quality of the point cloud, it is ensured that no two points are closer than r min h apart, and that every sphere of radius r max h in the domain has at least one point. Thus, the inter-point distance between each point and its nearest neighbor lies in the range (r min h, r max h). We follow conventionally used values of these parameters in Lagrangian meshfree GFDM literature, and set r min = 0.2 and r max = 0.4 [4,30]. This results in about 40 − 50 points in each interior neighborhood, with lesser at and near the boundary.
For the Lagrangian and ALE formulations, the movement of (parts of) the point cloud with the fluid velocity can result in the minimum and maximum inter-point distance criteria being violated. This happens in form of accumulation or scattering of points which would reduce the quality of the numerical results. To prevent this, points are added in holes containing insufficient points, and are merged in regions of accumulation. This method of fixing distortion is entirely local, and much cheaper than the remeshing done in mesh-based methods. Details about these procedures of adding and deleting points follow from [4,13,14,24,28].

Differential operators
GFDMs generalize classical finite differences to arbitrarily spaced point clouds, using a specialized weighted moving least squares approach. Consider a function φ defined on each point of the point cloud. At each point x j , numerical derivatives of φ are defined as a linear combination of function values in it's neighborhood where * = x, y, z, ∆, . . . denotes the derivative of interest, ∂ * is the continuous differential operator,∂ * j is the numerical differential operator at point x j , and φ l = φ (x l ). The numerical differential operators are thus given by the coefficients c * jl , which are independent of the function being differentiated. They are computed by a norm minimization process that ensures that monomials up to a specified order are differentiated exactly.
where M is the set of monomials being differentiated exactly. To compute the Laplacian, the monomials are complemented by the delta function to control the central stencil value c ∆ j j , which improves stability in the pressure Poisson equations [26]. In the present work, we consider monomials up to the order of 2. The weighting function W is defined such that neighboring points with the smallest distance to the considered point obtain the highest weight. In the present work, we use a truncated Gaussian weighting function for a constant c W > 0. We note that the same differential operators as defined above can also be equivalently derived by minimizing errors in Taylor expansions [26].
Using this procedure, we compute numerical gradient operators, and a numerical Laplacian. For more details on the computation of the differential operators, we refer to [4,14,26].

Time integration
3.3.1. Lagrangian formulation A strong form discretization of the physical model (Sect. 2) is done using the numerical differential operators defined above, and a chosen time integration scheme. For simplicity, the following considerations are based on a first order time integration.
Starting with the Lagrangian formulation, equations (1) can be rewritten as To improve readability, we henceforth use the shorthand ρ = ρ sol and c v = c v,sol . Together with equations (2)-(4), this is the starting point of the numerical discretization. The continuous spatial derivatives are replaced by their least squares approximated counterparts described in Sect. 3.2. We consider the superscript n + 1 to denote the next time level, and n for the current one, giving the time step size ∆t = t n+1 − t n . Below, we explain each of the steps of the discretization in the Lagrangian formulation for the microscopic scale simulations. Most of the steps are the same also for the macroscopic scale simulations, and the few differences are explained in Sect. 5.

Step 1. Point cloud movement
The discretization procedure begins by moving the point cloud according to a second order method [27] by with previous time step size ∆t 0 = t n − t n−1 .
Step 2. Temperature A semi-implicit time integration is then carried out to compute the new temperature T n+1 by with where the overhead ∼ indicates the discrete differential operators.
To simplify notation, the index of the points has been omitted. Equation (21) forms a sparse linear system of equations with unknowns T n+1 at each point of the point cloud. All sparse implicit linear systems arising in this and the coming steps are solved with a BiCGSTAB solver, without the use of a preconditioner.

Step 3. Concentrations
A similar procedure as that done for the temperature is carried out for the concentrations. We use a semi-implicit time integration for the concentration of each species c n+1 i , i = 1, . . . , N, with Step 4. ρ, η, and c v The updated density ρ n+1 , viscosity η n+1 sol , as well as heat capacity c n+1 v are then determined according to the definitions in Sect. 2.2.
Using the updated solution viscosity, a preliminary viscosity for the momentum equation is computed asη n+1 = η n+1 sol + η n turb .
Step 5. Hydrostatic pressure The pressure is split into its hydrostatic (body forces) and dynamic parts (movement of the fluid) as First the updated hydrostatic pressure p n+1 hyd is computed Using the updated hydrostatic pressure, a pressure guessp is computed which will be used while computing the new velocityp Step 6. Coupled velocity-pressure Time integration of the first equation in (19) provides the targeted divergence of velocity∇ T v n+1 . To solve for v n+1 and p n+1 in an implicit time integration scheme, we use the penalty formulation introduced in [10,13]. Using the pressure guess defined above, we obtain the following coupled velocity-pressure-system for preliminary velocityv n+1 and correction pressure p n+1 corr : and ∆t virt = A virt ·∆t, 0 ≤ A virt ≤ 1. If A virt = 1, the scheme corresponds to an implicit Chorin projection, see [3]. Theoretically, choosing A virt = 0 would give the exact solution. However, the linear system is ill-conditioned and can not be solved in most cases. For 0.001 ≤ A virt ≤ 0.1, conditioning of the linear system is sufficiently good. Furthermore, the resulting preliminary velocity features a divergence which is very close to the targeted one. We note that in equations (28) and (29), the stress tensor S n+1 was determined according to equation (2).
Step 7. Update velocity and pressure The updates of velocity and dynamic pressure are given by p n+1 dyn = p n dyn + p n+1 corr .
Step 8. Turbulence For the k-ε turbulence model, we derive a singularity formulation from equation (3): d dt If k, ε > 0 for all t n ≤ t ≤ t n+1 , numerical mean values can be determined from (31): We use the mean values to avoid singularities in the discretized k-ε turbulence model.
A fully implicit time integration scheme for the turbulent kinetic energy k n+1 can now be developed as A similar procedure is used to compute the updated turbulent dissipation The mean values are determined analytically. This is illustrated in detail for k ε m . Assuming that the diffusion term 1 ρ ·∆ η * k ε is negligible as well as defining we can rewrite equation (31) as Finally, the updated turbulent viscosity is determined by

Eulerian formulation
In case of the Eulerian formulation, [25] shows that a second order time integration scheme should be applied to numerically solve transport terms of the form v T ∇ in the GFDM context. For this purpose, the SDIRK2 method is proposed (see [1]), which features the same stability properties as an implicit Euler time integration scheme. Furthermore, an upwind discretization by means of a MUSCL reconstruction with a Superbee limiter is used.
The majority of the steps are the same as those carried out in the Lagrangian formulation. The movement step of the Lagrangian formulation is skipped here. And the coupled velocity-pressure system is modified to the following two-step procedure: 2 . Density and viscosity for the intermediate step can for instance be determined by linear interpolation between time levels n and n + 1.
In the second step, the preliminary velocity is determined aŝ

Further details
For more details on the Eulerian procedure, we refer to [25], and for similar GFDM Eulerian formulations, we refer to [29]. For numerical validations of the velocity-pressure scheme used here, their implementations within a GFDM framework, and a comparison of GFDM results with other numerical methods on benchmark problems, we refer to our earlier work [4,10,13,21,24,25,30].

MICROSCOPIC SCALE
To study the smaller scale (both spatially and temporally) dissolution of the salt species in the water, we consider representative geometries of the salt cavern in a so-called microscopic setup. In this section, we identify effective parameters of the dissolution process. Specifically, we compute the effective diffusion coefficient and the effective transition coefficient between water and surrounding species. These will be used later, in Sect. 5, in the macroscopic procedure to simulate the overall evolution of the salt cavern. The Lagrangian formulation is used here. The time integration of the underlying equations is done as presented in Sect. 3.3.1.
For the sake of brevity, we restrict the following description to sodium chloride as the species of interest. The same procedure can directly be transferred to any other species.

Setup
We consider a cylinder with diameter of 5m and height of 10m which is initially filled with pure water, i.e. c NaCl (t = 0) = 0. During the simulation, the temperature is fixed to T 0 = 20 • C.
The roof of the cylinder acts as an inexhaustible supply of sodium chloride which is modeled by applying a Dirichlet condition with saturation concentration For the hull of the cylinder, a homogeneous Neumann condition is applied. Aiming at a quasi-steady state, the bottom of the cylinder models an outflow boundary. In the interior, we solve where D micro = D laminar + D turb . The laminar diffusion coefficient for sodium chloride is given by D laminar = 1.611 · 10 −9 m 2 s (see [6]). For the turbulent part, we have Standard boundary conditions (Dirichlet and Neumann) are prescribed for velocity, pressure, and the turbulent quantities. The simulation runs until a quasi-steady state is reached, which will be explained below.

Evaluation strategy
In order to determine the effective quantities, the cylinder is split in the axial direction (z-direction) into equal sub-cylinders SC j , j = 1, . . . , J. These are used to estimate the mass flow. The planes between the sub-cylinders are denoted by help-planes HP j , j = 1, . . . , J − 1.
We note that the moving Lagrangian nature of the simulations means that point locations are always changing in each time step, except in the trivial case when v = 0 which does not occur here. Thus, a true steady state never occurs. Rather, simulations run till a quasi-steady state is reached, which is determined by the averaged values of the mass flow in the sub-cylinders. A quasi-steady state is said to be reached when the relative change of the mass flow in each of the sub-cylinders is within a tolerance specified (here, 10 −4 ) for 5 consecutive time steps.

Effective diffusion coefficient The mass flow of sodium chloride is given by
wherec NaCl is the mean concentration. The mass flow and the mean concentration in sub-cylinder SC j are determined by Based on the mean concentration in a sub-cylinder SC j , we can approximate its normal derivative with respect to the help plane HP j . This yields the effective diffusion coefficients in each sub-cylinder Once a quasi-steady state is reached, an overall effective diffusion coefficient can be determined. To accommodate the "quasi-steady" character of the simulation, we use a time-averaged effective diffusion coefficient, over a small time interval, and over each of the sub-cylinders. This value will later be used in the macroscopic setup.

Effective transition coefficient
The effective transition coefficient γ NaCl,eff is derived in a manner similar to that done for the effective diffusion coefficient D NaCl,eff above.
Once again, the time-averaged values of the effective transition coefficient in each of the sub-cylinders at the quasi-steady state gives the overall effective transition coefficient which will be used in the macroscopic simulations in Sect. 5.

Effective solution rate
With the help of γ NaCl,eff , we can define the solution rate of sodium chloride for given temperature T 0 by

Numerical results
In the simulations carried out, we choose J = 10 to divide the cylinder domain considered into 10 sub-cylinders of height 1m each. We consider several levels of resolution to study the convergence of the effective parameters being determined to resolution-independent values. The coarsest resolution used is h = 0.8m corresponding to 70400 points in the domain. h is consecutively halved till h = 0.1m corresponding to 20892871 points in the domain. Several resolutions in between are also considered to better illustrate the converged values of the effective parameters. We note that the number of points mentioned here are at the initial time of the simulation (t = 0). This number of points will slightly vary in time due to the addition and deletion of points explained in Sect. 3.1. The evolution of the concentration for h = 0.18m is illustrated in Fig. 2 in the time interval [0s, 100s]. As expected, the flow is characterized by viscous fingering.   The convergence of effective diffusion as well as transition coefficient with decreasing h is shown in Fig. 3 and Fig. 4, respectively. The values plotted are also tabulated in Table II, along with the relation between the interaction radius h and the number of points in the domain NP. The time step size is governed by ∆t = CFL Lag · h |v| , with CFL Lag set to 0.2. The diffusion coefficient converges to D NaCl,eff = 0.1 m 2 s , while the transition coefficient converges to γ NaCl,eff = 0.000042 m s . Using equation (53), we obtain a maximum solution rate of R NaCl,max (20 • C) = 0.0150 kg m 2 ·s . Compared to the solution rate of 0.0488 kg m 2 ·s for T 0 = 23 • C determined in [11] at a crystal level, the estimated solution rate is of the correct order of magnitude.

MACROSCOPIC SCALE
We now model the overall evolution of the salt cavern during the double-well solution mining process. Both the Lagrangian as well as the Eulerian formulation are evaluated for this.
The model equations and time integration procedures for the macroscopic scale simulations are the same as those described in Sect. 2 and Sect. 3.3, respectively, with a few variations. Firstly, the dissolution of the salt into the water occurs at much smaller spatial and temporal scales than those used here. To take this into account, the dissolution process of the salt at the cavern walls are modeled using a Robin boundary condition for the concentration Here, the effective diffusion coefficient D NaCl,eff , as well as the effective transition coefficient γ NaCl,eff are the values determined in the microscopic simulation in Sect. 4, D NaCl,eff = 0.1 m 2 s and γ NaCl,eff = 0.000042 m s . A further difference in the time integration procedure comes in Steps 2 and 4 described in Sect. 3.3. Here, we fix the temperature to T 0 = 20 • C and, subsequently, obtain the corresponding saturation concentration c s NaCl = 357 kg m 3 . For simplicity, the following linearized relations for density and viscosity of the solution are used (see [25]) η(c NaCl ) ≈ (1.96 · 10 −6 · c NaCl + 10 −3 ) Pa s .

Setup
We are interested in the geometrical evolution of the double-well salt cavern. The initial geometry is given by a small cavern filled with pure water that is surrounded by sodium chloride, see Fig. 5. The dimensions of the initial cavern are approximately: width of 90m, height of 50m, and depth of 26m. The sodium chloride deposit is limited to impermeable surrounding rock. The pipe on the left side acts as an inlet of fresh water with inflow velocity |v in | = 1 m s , whereas the pipe on the right side acts as the outlet.
In reality, the maximum diameter of the pipes is of the order of 1m. Hence, the resolution of the point cloud close to the inlet and the outlet has to be of the order of 0.1m to ensure accurate results in case of the Lagrangian formulation. This would lead to an extremely small time step size compared to the desired simulation time of several years/decades due to the CFL-condition Numerically, we observe that stable results are achieved for CFL Lag = 0.15, which results in ∆t Lag = O(0.1s). Performing long-term simulations over months of simulation time is not feasible with such a small time step. This results in the need for using the Eulerian formulation for the problem at hand.
In order to allow for a comparison of Lagrangian and Eulerian formulation, we consider pipes of diameter 12m. This also decreases the required actual time being simulated. Numerically, we observe Figure 5. Macroscopic simulation setup -initial geometry, see [25]. that with the significantly larger diameter used here, the evolution of the cavern only requires a few hours of physical time to be simulated, compared to the few months or years with the actual diameter. We note that despite this time reduction, this still corresponds to a time scale two orders of magnitude greater than that used in the microscopic simulations in Sect. 4.
The simulations are performed with a constant interaction radius of h = 4m. An important point to note here is that this spatial resolution considered is of the same order of magnitude as the height of the sub-cylinders in the microscopic simulations in Sect. 4.

Movement of the boundary
The movement of the boundary of the cavern can be defined by the Stefan condition ρv = γ NaCl,eff (c s NaCl − c NaCl ), see [9]. This yields and, consequently, a movement of the boundary in normal direction n with velocity v boundary = v · n.
To speed up computation, a time lapse procedure can be applied [25]. Due to small flow velocities inside the salt cavern, an additional speed-up factor A can be introduced in the definition of v by For stability reasons, we require The maximum movement velocity of the boundary occurs in case of pure water, i.e. c NaCl = 0, and is given by Given a maximum time step size ∆t max , equation (60) leads to the constraint Figure 6. Translational velocity in the macroscopic simulation, see [25] (Eulerian formulation including ALE at the moving boundary).
Due to the movement of the boundary, interior points close to this boundary have to move in the Eulerian formulation also. For this purpose, the ALE-approach presented in [8] is used. Based on current and future position of an affected interior point, the translational velocity is determined, see Fig. 6. Due to the explicit movement of these points, the convection terms in the numerical model in Eulerian form in Sect. 3.3 must refer to the relative velocity v − v trans instead of v. Furthermore, this introduces a CFL-condition of the form This depends on the boundary velocity v = O(0.01 m s ) which is considerably smaller than the flow velocity (the inflow velocity is 1 m s , while the maximum velocity in the domain is even bigger). h min is subject to the desired resolution at the moving boundary. At the inlet and the outlet, a coarse resolution is sufficient in this case.

Numerical results
Starting from the initial domain as shown in Fig. 5, the salt cavern expands till the outer domain is filled. Physically, this outer domain can represent either a rock formation where the salt cavern ends, or prescribed limits of the region where the solution mining is to be carried out. Fig. 7 illustrates the evolution of the salt cavern in the Eulerian formulation according to C = c NaCl . An animation of this process can be found in the Online Resource. This expansion is quantified by plotting the volume as a function of time in Fig. 8. We note that once the entire outer domain is filled, the volume of the domain about 27.4 times that of the initial domain. This shows the need of a meshfree method for the present application. If a mesh-based method were to be used for this simulation, the entire domain would need to be meshed initially, and an expensive and less accurate tracking of the expansion would need to be carried out.
To compare the results of the Eulerian and Lagrangian formulations, we consider simulations on the same initial point cloud with h = 4m which corresponds to NP 0 = 90 744 points at the initial state. ∂ Ω out c NaCl v · n dA dτ , where ∂ Ω out is the outflow boundary located at the top of the extraction well. Physically, this represents a measure of the concentration of salt being extracted. The time evolution of Q c is shown in Fig. 9. It illustrates that both formulations produce very similar results.
To emphasize the need of the ALE formulation for such a simulation, we compare the time steps required in both the ALE and Lagrangian formulations for stability. Considering the simulation time of 7200s, we observe that the Lagrangian formulation required at least 22915 time steps to obtain stable results, which corresponds to an average time step size of ∆t ≈ 0.31s. On the other hand, similar results, as shown in Fig. 9, can be obtained in the Eulerian formulation (with ALE near the boundaries) with only 936 time steps corresponding to an average time step size of ∆t ≈ 7.69s, which is approximately 25 times that needed in the Lagrange case. on a macroscopic, as well as a microscopic scale. Both Lagrangian and Eulerian approaches were considered.
On the macroscopic scale, we considered the expansion of the salt cavern as a result of erosion occurring as the salt dissolves in the water. In reality, this procedure occurs over the time span of several months or years. A simplified geometry was considered here, which enabled a comparison between the Eulerian and Lagrangian formulations. In this simplified macroscopic set-up, the expansion of the salt cavern occurred over the time scale of several hours. Since the dissolution of salt in water occurs on a much smaller time level we also considered a microscopic set-up over a duration of a few minutes. This was used to determine effective parameters governing the dissolution process. Using the example of sodium chloride as the species of interest, effective diffusion and transition coefficients were determined in the microscopic simulations. These values were then used in the macroscopic simulations to determine the evolution of the concentration inside the salt cavern and to specify the solution rate of the salt species at the boundary, i.e. to model the geometrical evolution of the salt cavern.
A comparison of the numerical results of the Lagrangian and Eulerian formulations (extended by an ALE-approach) in the macroscopic case illustrates the advantages of the latter one due the possibility of using much larger time step sizes. Aiming at a simulation time of several years, the forecast computation time for a simulation of a double-well solution mining process based on the Lagrangian formulation would be of the order of years. In contrast to that, the flexibility of the Eulerian formulation regarding the resolution of the point cloud (local refinement only at the moving boundary) enables meshfree simulations in reasonable time -especially in terms of real applications.

CONFLICT OF INTEREST
On behalf of all authors, the corresponding author states that there is no conflict of interest.