Geomechanical modelling of sinkhole development using distinct elements: model verification for a single void space and application to the Dead Sea area

Abstract

Mechanical and/or chemical removal of material from the subsurface may generate large subsurface cavities, the destabilisation of which can lead to ground collapse and the formation of sinkholes. Numerical simulation of the interaction of cavity growth, host material deformation and overburden collapse is desirable to better understand the sinkhole hazard but is a challenging task due to the involved high strains and material discontinuities. Here, we present 2-D distinct element method numerical simulations of cavity growth and sinkhole development. Firstly, we simulate cavity formation by quasi-static, stepwise removal of material in a single growing zone of an arbitrary geometry and depth. We benchmark this approach against analytical and boundary element method models of a deep void space in a linear elastic material. Secondly, we explore the effects of properties of different uniform materials on cavity stability and sinkhole development. We perform simulated biaxial tests to calibrate macroscopic geotechnical parameters of three model materials representative of those in which sinkholes develop at the Dead Sea shoreline: mud, alluvium and salt. We show that weak materials do not support large cavities, leading to gradual sagging or suffusion-style subsidence. Strong materials support quasi-stable to stable cavities, the overburdens of which may fail suddenly in a caprock or bedrock collapse style. Thirdly, we examine the consequences of layered arrangements of weak and strong materials. We find that these are more susceptible to sinkhole collapse than uniform materials not only due to a lower integrated strength of the overburden but also due to an inhibition of stabilising stress arching. Finally, we compare our model sinkhole geometries to observations at the Ghor Al-Haditha sinkhole site in Jordan. Sinkhole depth $/$ diameter ratios of 0.15 in mud, 0.37 in alluvium and 0.33 in salt are reproduced successfully in the calibrated model materials. The model results suggest that the observed distribution of sinkhole depth $/$ diameter values in each material type may partly reflect sinkhole growth trends.

Introduction

Sinkholes are enclosed surface depressions in sediments and rocks. They commonly result from subsidence of overburden into void space that is generated through the physical–chemical removal of material in the underground. In the final stage of a sinkhole process, a sudden collapse of the overburden may occur (Gutiérrez et al., 2014; Waltham et al., 2005). Removal of material and void formation in the underground is usually related to hydraulic flow and to associated dissolution, physical erosion of material or both. Subsidence may occur continuously over a long time depending on the flow conditions and material properties (Goldscheider and Drew, 2007; Parise and Gunn, 2007; Waltham et al., 2005). Depending on the properties of the overburden (cover or caprock) and the evolution stages, different sinkhole morphologies can be described. Typical end-members can be defined (Fig. 1; see Gutiérrez et al., 2008, 2014).

Conceptual models of sinkhole formation. (a) Subsurface dissolution and/or subrosion caused by focussed flow in soluble/weak material. (b) Cover suffusion or dropout sinkhole that forms by material transport through a pipe or along a funnel. A weak cover material slumps into the voids and creates a sinkhole with low depth $/$ diameter (De $/$ Di) ratio and flat-to-steep margins depending on the material cohesion. (c) Cover or caprock collapse sinkhole. Large voids may stay initially stable in a strong material, but their growth leads to a sudden overburden collapse. The formed sinkholes have usually a high De $/$ Di ratio and contain steep margins with large ground cracks. Both sinkhole types represent late-stage end-members and mixtures of both are very common in nature (Fig. 2).

[Figure omitted. See PDF]

Sinkhole development at the Dead Sea

The Dead Sea is a hypersaline terminal lake and is one of the world's most active areas of sinkhole development. More than 6000 sinkholes have formed there at an increasing rate over the last 35 years (Abelson et al., 2017). Previous studies relate the sinkhole formation at the Dead Sea to the regression of the lake, which has been ongoing since the 1960s, and to the consequent invasion of evaporite-rich sedimentary deposits around the Dead Sea by relatively fresh groundwater. Evaporitic minerals in the sediments are susceptible to dissolution, while the non-evaporitic sedimentary materials are weak (poorly consolidated or unconsolidated) and can easily be physically eroded by subsurface flow (“piping”). Some studies have highlighted the role of subrosion, i.e. both mechanical and chemical (leaching) erosion of the subsurface (Wadas et al., 2016), in the development of sinkholes (Al-Halbouni et al., 2017; Arkin and Gilat, 2000; Polom et al., 2018), while others have focussed on the role of dissolution only in generating large cavity development in a relatively shallow but thick salt layer (Ezersky and Frumkin, 2013; Taqieddin et al., 2000; Yechieli et al., 2006).

In this paper, we draw upon observations from the sinkhole site of Ghor Al-Haditha (31${}^{\circ }$18${}^{\prime }$45${}^{\prime \prime }$ N, 35${}^{\circ }$31${}^{\prime }$52${}^{\prime \prime }$ E) on the southeastern shore of the Dead Sea in Jordan (Figs. 2, 3). Sinkhole development in the area has been active since 1986 (Sawarieh et al., 2000), with ongoing damage or destruction of infrastructure and agriculture. As of 2018, the cumulative number of sinkholes formed there has passed 1000 (Watson et al., 2018). Photogrammetric datasets have been acquired in 2014, 2015 and 2016 to produce high-resolution and high-accuracy digital surface models and orthophotos for the Ghor Al-Haditha sinkhole site. Although the results for 2015 and 2016 shown below in this paper are new, the methodology of their generation is the same as for the 2014 survey, which is described in detail by Al-Halbouni et al. (2017).

Sinkhole examples from the eastern shoreline of the Dead Sea. (a) Sinkhole with De $/$ Di $\sim \mathrm{0.15}$ formed in semi-consolidated lime-carbonate mud of the former Dead Sea bed. (b) Sinkhole with De $/$ Di $\sim \mathrm{0.33}$, formed in semi-consolidated sandy gravel (alluvial) sediment. Note the deep cracks and tilted blocks marginal to the sinkhole. (c) Sinkhole with De $/$ Di $\sim \mathrm{0.33}$ formed in semi- to well-consolidated salt material. Note the typical overhanging sides and pronounced cracks in the surroundings.

[Figure omitted. See PDF]

Sinkholes form in the three “end-member” near-surface materials at the Ghor Al-Haditha sinkhole site (Fig. 2): (1) unconsolidated to semi-consolidated lacustrine clayey carbonates (“mud”) with interleaved thin evaporite layers; (2) unconsolidated to semi-consolidated alluvial sand–gravel sediments; and (3) rock salt (mainly halite) with interleaved thin mud layers. The main morphological distinction is that narrower and deeper sinkholes occur in the “alluvium” and in the “salt” (Fig. 3b, c), whereas wider and shallower sinkholes occur in the “mud” (Fig. 3a). Many sinkholes in the alluvium and especially in the salt have overhanging sides and/or large marginal blocks and deep (up to several metres) concentric ground cracks. The alluvium and the salt can sustain metre-scale or multi-metre cavities associated with sinkhole development (Al-Halbouni et al., 2017; Closson and Abou Karaki, 2009; Yechieli et al., 2006). The mud sinkholes commonly contain a wide peripheral zone of back-rotated blocks delimited by small faults that down-throw towards the centre. Ground cracks are commonly also well developed around sinkholes in the mud but are not as deep (up to a few tens of centimetres) as in the other materials.

Representative morphological data from single sinkholes at the eastern shoreline of the Dead Sea: (a) in the mud flat, (b) in the alluvial material and (c) in the salt cover. Shown here are orthophotos (left column), digital surface models (DSMs, middle column) and topographic cross-sections (right column) with a resolution of 10 cm and an accuracy of 12 and 17 cm ($H$ and $V$), respectively. These were created from low-altitude aerial images acquired in 2015 and processed by structure-from-motion photogrammetry. Contours of elevation in metres are indicated for clarity on the DSM, which is plotted in the same colour scale for all materials.

[Figure omitted. See PDF]

Numerical modelling of sinkhole development

The numerical simulation of sinkhole development is of interest to understand better the processes of sinkhole formation and the related hazard. Continuum-mechanics approaches (Carranza-Torres et al., 2016; Fazio et al., 2017; Fuenkajorn and Archeeploha, 2010; Parise and Lollino, 2011; Rawal et al., 2016; Salmi et al., 2017) have generally defined a single cavity in an elastic or elastoplastic half space and assessed the static threshold strength of the overburden to predict mechanical failure. This approach is possibly suitable for assessing the factor of safety of an individual, fully developed cave and for deriving a relation between measured surface subsidence and cavern configuration. However, the geometries of voids involved in sinkhole development are often non-singular, irregular and distributed on lots of scales (Abelson et al., 2017; Al-Halbouni et al., 2017; Ezersky et al., 2017; Gutiérrez et al., 2016; Parise et al., 2018; Yizhaq et al., 2017). Alternatively, continuum-based corrosion models have addressed the rock dissolution and void growth in a hydrogeological framework (Kaufmann and Romanov, 2016; Shalev and Lyakhovsky, 2012). This approach has the advantage of accounting for geometrically complex or stochastic void development and the role(s) of material heterogeneity, but it does not account for effects of overburden instability. Both of these past continuum-based approaches have neglected the mechanical consequences of void growth and the explicit simulation of sinkhole collapse.

Distinct element method (DEM) modelling is increasingly used in geoscience for numerical simulation of large-strain and discontinuous rock deformation (Cundall and Strack, 1979; Potyondy and Cundall, 2004). The main advantage of the DEM is its ability to simulate rock samples or rock masses as an assemblage of discrete particles or blocks, which can undergo large displacements and rotations. The method uses a so-called soft contact approach where the particles are rigid but allowed to overlap at contact points. Based on updated particle positions, the contacts between particles are automatically detected or deleted during the simulation. Based on the relative displacement and velocity of the particles in contact, interaction laws are used to update the forces and moments transmitted through the contacts. The resultant forces and moments that accumulated on each particle are subsequently used to solve Newton's second law of motion and to update the particle's position and velocity. Elastoplastic bonds of finite strength can be accounted for in the interaction law and enable a quasi-continuum behaviour at assembly scale, which can evolve to highly discontinuous deformation as bonds between particles break and damage develops. In this way, the DEM can overcome limitations of continuum-based numerical simulation of large and highly localised strains in discontinuous media (Jing and Stephansson, 2007). Using the DEM, recent advances have been made in, for example, rock mechanics (Schöpfer et al., 2009), slope stability and mass movements (Thompson et al., 2010), mine or tunnel stability (Bonilla-Sierra et al., 2012), synthetic rock mass modelling (Ivars et al., 2011), fracture growth (Schöpfer et al., 2016), hydrofracture and caldera subsidence analysis (Holohan et al., 2011, 2015, 2017).

For modelling sinkhole collapse, Baryakh et al. (2008, 2009) used the DEM to conduct simple stability tests and mechanical analyses for a single, instantaneously generated cavity of varying geometry, depth and overburden mechanical properties. Other studies have adopted a similar approach but also included discrete fracture networks (DFNs) that represent predefined or empirically determined discontinuities (joints, faults) within rock masses (Hatzor et al., 2010). DEM coupled with finite element modelling (FEM) has been used for simulating mechanical failure above a large salt cavity (Mercerat, 2007), but the DEM part was limited to a single rock layer within the overburden. Again, the main shortcoming of these earlier DEM-based studies is that cavity growth and related mechanical development were not explicitly simulated.

Contribution of this paper

This paper reports the first two-dimensional DEM simulations of sinkhole formation that explicitly simulate both void growth and overburden collapse. In part, the approach of void growth adopted here is similar to that in a recent work on mine caving (Sainsbury, 2012). Our study builds upon the previous works of Caudron et al. (2006) and Baryakh et al. (2008, 2009) but goes further in calibration of geomechanical behaviour, in complexity of process and in application to natural sinkholes. As in many previous studies, we focus here on the creation, growth and instability of a single void leading to a single sinkhole. We use the general-purpose commercial DEM Particle Flow Code in two dimensions (PFC2D) software, developed by Itasca Consulting Group Inc. (Potyondy, 2014b). Further details of the DEM, as implemented in PFC2D, are covered in Appendix A. Note that, in accordance with PFC2D convention, compressive stress is taken as negative throughout this paper.

Regarding the structure of this paper, we begin by summarising tests on model resolution, model dimensions and void creation procedures, and we show results of benchmarking to continuum-based solutions for displacement around a void. We then show the results of calibration tests that were used to tune the bulk geomechanical behaviours of the DEM particle assemblies. Following this, we analyse the evolution stages of model void growth and sinkhole collapse for uniform and layered materials. We then compare the morphological parameters at the Ghor Al-Haditha survey site to those predicted by our models. In the final part, an outlook to future improvements and applications is given.

A distinct element method approach for modelling cavity and sinkhole formation

In this section, we report on convergence and benchmarking tests for the DEM model as pertained to cavity generation. To this end, we firstly simulate a material that behaves elastically by using bond strength and cohesion values at the upper limit of realistic rock strengths (see Table 1). We also report on the material parameter calibration by simulated biaxial compression and tension tests applied to the numerical materials mimicking those common in the Dead Sea region. Finally, we summarise the final procedure for cavity growth that is based on these tests but implemented under conditions in which the DEM model material is weak enough to fail and lead to sinkhole formation.

Setups for model benchmarking, calibration and sinkhole simulation. Panel (a) indicates model verification and benchmarking. A circular cavity of radius $r=\mathrm{5}$ m is placed instantaneously at a depth $h=\mathrm{35}$ m and a distance $d=H-h$ from a fixed point at the bottom of a box of varying dimensions $H\times W$ and particles of radii $R$ subject to the body force due gravitational acceleration $F_{\mathrm{g}}$. Panel (b) indicates rock tests for material parameter calibration. A sample is contained within walls that are used for applying confining pressure simulating the materials' response at different depth. A servo-mechanism controls the walls' axial velocity. For a tensile test, grips of certain thickness are defined at the bottom and top of the sample and moved outwards. (c) Sinkhole simulation by quasi-static incremental single void growth. $T/D$ is referring to the overburden thickness to diameter ratio of either a stable cavity ($D=D_{\mathrm{cav}}$) or an unstable collapse zone ($D=D_{\mathrm{col}}$). Yellow/red circles represent particles that act as extensometers/markers, respectively. Big red circles indicate overlapping measurement circles distributed within an area of interest.

[Figure omitted. See PDF]

Determination of the optimal void space installation and model dimensions

We tested model sensitivity to resolution, dimension and void installation method. In a first test, different void space installation methods were compared in terms of computation time. For this, surface particle displacement was tracked above a cavity of 5 m radius placed at 35 m depth (Fig. 4a), an assumed realistic subrosion zone depth and dimension. Two methods utilised a particle deletion scheme, while two other methods were based on particle radii reduction. No substantial differences in the vertical and horizontal surface displacements were observable; i.e. the methods did not affect the outcome of the elastic solution, but the particle deletion scheme was 1 order of magnitude faster than radii reduction. Hence, the particle deletion scheme was chosen as appropriate for the following tests and the sinkhole models. More details on the results and the set of investigated parameters can be found in Appendix B1.

In a second test, model width, height and particle radii were varied to determine the optimal model dimensions for the problem of a void space in the subsurface. The void installation method based on instantaneous particle deletion was applied. The final results indicates that symmetric boundaries of $H\times W=\mathrm{400}\times \mathrm{400}$ m with a particle mean radius of 0.32 m yield the best results. These model dimensions and resolution were hence chosen for the main model set reported below. Details on the convergence tests that led to this choice of dimensions and resolution can be found in Appendix B2.

Benchmarking of the DEM approach against analytical solutions and BEM

We performed a benchmarking of surface displacements in the DEM cavity development models with displacements derived from different continuum-based approaches. Cavity depth and size, model dimensions (Fig. 4a) and the bulk elastic parameters of the DEM material in Table 1 serve as input parameters for two analytical solutions and a boundary element (BEM) numerical model.

Bulk properties of the particle assemblies used in the benchmarking of DEM cavity formation models vs. analytical solutions and BEM.

Bulk parameter (unit)	Symbol	Unit	Value
Porosity	$n$	–	0.16
Density${}^{\mathrm{1}}$	${\mathit{\rho }}_{\mathrm{bulk}}$	(kg m${}^{-\mathrm{3}}$)	2100
Earth pressure coefficient at rest${}^{\mathrm{1}}$	$K_{\mathrm{0}}$	–	0.26
Young's modulus	$E$	(GPa)	$\mathrm{5.337}\pm {\mathrm{7}}^{-\mathrm{3}}$
Poisson ratio${}^{\mathrm{2}}$	$\mathit{\nu }$	–	$\mathrm{0.39}\pm \mathrm{0.03}$

${}^{\mathrm{1}}$ Estimation based on measurement circles in five different particle assemblies at a depth of 35 m. ${}^{\mathrm{2}}$ Estimation based on fitting of measurement circle data in 10 different simulated rock tests. Compare Sect. 2.4.

The analytical solutions are for a circular cavity in a gravity loaded, infinite, linear elastic full/half space under plane strain conditions (Kirsch, 1898; Verruijt and Booker, 2009). The Kirsch solution is a classical full-space solution for simple excavation shapes but does not include the free-surface effect; mathematical details are provided by Brady and Brown (2006). The Mindlin (1940) solution is for stresses around tunnels and includes free-surface effects; mathematical details are given by Verruijt and Booker (2009). The input values for the Mindlin analytical solution are $d/h=\mathrm{4}$, $E=\mathrm{5}$ GPa, $\mathit{\nu }=\mathrm{0.39}$ and $K_{\mathrm{0}}=\mathrm{0.26}$. The reader is referred to Appendix B3 for more details on the effect of $d$ and $E$.

The BEM model is based on a code by Nikkhoo and Walter (2015) and simulates the surface displacements along a cross-section above a 3-D cylindrical void space. The void space is simulated as a traction-free, horizontal, north–south-oriented cylinder of 200 m total length. The cylinder's centroid is located exactly beneath the origin; a hydrostatic remote stress is applied equal to the gravitational stress ${\mathit{\sigma }}_{xx}={\mathit{\sigma }}_{yy}={\mathit{\sigma }}_{zz}=\mathit{\rho }gh$, where $h$ is the depth to the cylinder centroid.

The DEM model displacements $\left(U_x,U_y\right)$ as well as the displacement differences $\left(\mathrm{\Delta }{U}_x,\mathrm{\Delta }{U}_y\right)$ in Fig. 5 match closely the Kirsch solution and the BEM results. For the Mindlin solution, this is only true for the horizontal components. For the vertical components, the modelled components only match the Mindlin solution in the near field of the subsidence centre, while in the far field a large disagreement is observable, expected from an intrinsic mathematical difficulty in determining the displacement of a stress-loaded half space (Appendix B3).

Results of benchmarking of the DEM cavity model against continuum-based cavity models: analytical solutions and modelled displacement curves for model dimensions of $\mathrm{400}\times \mathrm{400}$ m with particle radius 0.32 m. The reference is a cavity with $r/h=\mathrm{0.143}$, at 35 m depth and with a 5 m radius. (a) Horizontal displacement, (b) vertical displacement, (c) horizontal displacement difference and (d) vertical displacement difference.

[Figure omitted. See PDF]

Model setup for cavity generation and sinkhole formation

Based on the above-described tests for model resolution, dimensions and cavity generation, a generalised setup for cavity growth with attendant fracturing and sinkhole collapse is presented in Fig. 4c. This setup comprised a $\mathrm{400}\times \mathrm{400}$ m assembly of parallel-bonded particles of 0.32 m radius on average. The assembly is subdivided according to bond and particle contact properties into a cover material sequence that lies over a “soluble” or “mobile” material with a fixed basement rock below. The cavity is grown according to a material removal zone of arbitrary geometry, taken here as a vertically orientated half ellipse, using the incremental deletion approach (M2) described above. For the technical details of the bond installation and void space creation procedure, see Appendix B5.

A key geometric parameter in subsidence studies is the ratio of overburden thickness ($T$) to width or diameter ($D$) of the undermined area before the initiation of subsidence or collapse. In the following, $T/D$ refers to (1) thickness of the overburden/diameter of the cavity ($T/D_{\mathrm{cav}}$) if materials can sustain a cavity in or around the material removal zone or (2) thickness of the strong material/diameter of the destabilised zone ($T/D_{\mathrm{col}}$) if materials cannot sustain a cavity. This subsurface destabilised zone is shown arbitrarily as a triangular-shaped area in Fig. 4c. For each model setup, at least two (mostly five) different particle assemblies were run, and the errors in the following are based on these.

Calibration of the distinct element method approach for modelling sinkhole collapse at the Dead Sea

Bulk parameter estimation of geomaterials adjacent to the Dead Sea

The geotechnical parameters of rocks and soils commonly cover a wide range, as they depend strongly on detailed mineral composition, grain sizes, external stress conditions, fluid saturation and stress histories; cf., e.g. Brady and Brown (2006) and Jaeger et al. (2007). Here, we consider geotechnical parameters for the three main material types involved in sinkhole formation at the Dead Sea region: (1) lacustrine clayey carbonates and evaporites, (2) alluvial sands and gravels and (3) pure rock salt (halite) (Table 2).

Estimated geomechanical properties of main materials in sinkhole-affected areas at the Dead Sea. References: Ezersky et al. (2017); Ezersky and Livne, (2013); Frydman et al. (2008, 2014); Hoek (2007); Khoury (2002); Manger (1963); Polom et al. (2018); Zhu (2010).

Parameter	Symbol	Unit	Wet lacustrine mud	Alluvial sediment	Holocene salt
Bulk density	${\mathit{\rho }}_{\mathrm{bulk}}$	(kg m${}^{-\mathrm{3}}$)	1500–2100	1500–2300	1400–2200
Porosity	$n$		0.2–0.3	0.2–0.3	0.05–0.4
Friction angle	$\mathit{\varphi }$	(${}^{\circ }$)	2.4 (wet)–34 (saturated)	30–40	53
Young's modulus	$E$	(MPa)	83	220	300–10 000
Poisson's ratio	$\mathit{\nu }$		0.2–0.4	0.15–0.35	0.2–0.4
Unconfined compressive strength	UCS	(MPa)	$-\mathrm{0.243}$ to $-\mathrm{0.053}$	$-\mathrm{5}$ to $-\mathrm{0.1}$	$-\mathrm{5}$ to $-\mathrm{1}$
Unconfined tensile strength${}^{*}$	$T$	(MPa)	0.0053–0.0243	0.01–0.5	0.1–0.5
Cohesion	$c$	(MPa)	0–0.019	0.027–1.33	0.84–1.73
					(at 20–40 m depth)

${}^{*}(T\sim$ UCS/10) after Hoek (1968).

For lacustrine mud, friction angle, cohesion, porosity and density parameters from laboratory tests are used (Ezersky et al., 2013, 2017; Frydman et al., 2008, 2014). For the alluvial sediments, upper limits are given by nearby field investigations in firm sandstone rocks (El-Naqa, 2001) and also by published values for medium-grained Quaternary sand–gravel (Carter, 1983; Manger, 1963; Taqieddin et al., 2000). The bulk modulus of alluvial sand–gravel and lacustrine clays were estimated using Poisson's ratio values from the literature (Zhu, 2010) and shear-wave velocities from recent field measurements (Polom et al., 2018), where the latter were reduced by a factor of 1.5 to account for drained conditions.

Elastic parameters and strength values of the field materials have been estimated by using tables from Brown (1981) and Hoek (2007) and by classifying the clayey mud as grade R0 in terms of intact rock consistency and the alluvial sediments as grade R0–R1. The Holocene salt rock of the Dead Sea is considered weaker than typical halite rock salt (Frydman et al., 2008, 2014) and has been classified as grade R1. The cohesion value of salt is strongly depth dependent and has been determined by using depth-normalised results derived from triaxial tests (Frydman et al., 2014) via $$c={\displaystyle \frac{q\times z}{\mathrm{2}\times N_f}},$$ where $N_f=\sqrt{\frac{\mathrm{1}+\sin \mathit{\varphi }}{\mathrm{1}-\sin \mathit{\varphi }}}$, $c$ the cohesion, $z$ the depth of the rock sample, $q$ the intercept in a principal stress ${\mathit{\sigma }}_{\mathrm{1}}\left({\mathit{\sigma }}_{\mathrm{3}}\right)$ plot and $\mathit{\varphi }$ the friction angle. We use a friction angle of $\mathit{\varphi }=\mathrm{54}$${}^{\circ }$, depth $z=\mathrm{20}$–40 m and an intercept of $n=\mathrm{259}$ kPa for a specific rock weight of 18 kN m${}^{-\mathrm{3}}$ (Frydman et al., 2014). For the alluvial sediment, we assume a friction angle of $\mathit{\varphi }=\mathrm{34}$${}^{\circ }$ and an UCS of 0.1–5 MPa and calculate the cohesion value by the well-known relation $c=\mathrm{UCS}/\left(\mathrm{2}{N}_f\right)$ (Jaeger et al., 2007). Modulus, friction angle and strength hereby depend strongly on the porosity distribution, while Poisson's ratio is quasi-independent of it (Schöpfer et al., 2009).

Calibration of DEM material properties via simulated rock tests

Bulk material parameters are determined by simulated, biaxial compression and tension tests similar to those described by Khanal and Schubert (2005) and Schöpfer et al. (2007) (see Fig. 4b). We generated material “samples” with dimensions of $\mathrm{10}\times \mathrm{8}$ m and with a mean particle radius of 0.32 m and an initial porosity of 0.2. Each sample then contains approximately 200 particles. In order to simulate the materials of the Dead Sea region, we used the microproperties listed in Table 3 and the assembly generation scheme outlined in Appendix B5.1. Note that, for the lacustrine mud material, we also implemented a re-bonding (“annealing”) scheme to simulate the cohesiveness of that material. Tests were conducted with confining pressures $p$ between 0 and $-\mathrm{5}$ MPa, corresponding for a bulk density of $\sim \mathrm{2000}$–2200 kg m${}^{-\mathrm{3}}$ to depth range of 0–250 m. Measurement circles (averaging regions as described in Potyondy and Cundall, 2004) are installed in the centre of the sample to determine the stress–strain and porosity values.

Particle, contact and bond properties for DEM sinkhole collapse models. Note the geology convention.

Parameter${}^{\mathrm{1}}$	Symbol	Unit	Lacustrine mud	Alluvial sediment	Holocene salt
Initial material porosity	$n$	–	0.2	0.2	0.2
Particle density	$\mathit{\rho }$	(kg m${}^{-\mathrm{3}}$)	2715	2750	2500
Contact Young's modulus${}^{\mathrm{2}}$	$E_L$	(GPa)	0.1	0.2	1
Contact Young's modulus at	$E_{\mathrm{W}}$	(GPa)	5	5	5
particle-wall contacts
Bond Young's modulus	$E_{\mathrm{Pb}}$	(GPa)	0.1	0.2	1
Bond tensile strength	${\mathit{\sigma }}_c$	(MPa]	0.1	0.5	1.0
Bond cohesion	ć	(MPa)	0.02	0.5	1.0
Bond friction angle	$\mathit{\varphi }$	(${}^{\circ }$)	2.4	34	54

${}^{\mathrm{1}}$ See Appendix A and Potyondy and Cundall (2004) for a detailed definition of these parameters. ${}^{\mathrm{2}}$ A friction coefficient of 0.5 and a normal-to-shear stiffness ratio of 2.5 is chosen for all materials.

The sandy gravel and salt materials show brittle failure behaviour (i.e. a sharp post-peak stress drop) at low confining pressures, which changes to brittle–ductile behaviour for larger confining pressures (Fig. 6). Ductile is defined as the state of deformation without significant loss of strength, and the transition to this behaviour is the brittle–ductile transition (Byerlee, 1968). The salty mud material shows a brittle–ductile transition for all tested confining pressures or, more precisely, a brittle-to-cataclastic-flow transition to distinguish it from the brittle-to-crystal-plastic transition (Schöpfer et al., 2013).

Differential stress vs. vertical strain plots for simulated compression tests. Confining pressures of $-\mathrm{0.1}$, $-\mathrm{1.5}$ and $-\mathrm{4.0}$ MPa are compared that highlight the depth-dependent division into elastic, yielding and post-peak behaviour in all tested materials: (a) low-strength lacustrine mud, (b) middle-strength alluvial sediment and (c) high-strength Holocene salt material. Black dashed line in panel (a) marks the cutoff limit at $E_{yy}=-\mathrm{0.02}$ for lacustrine mud elastic properties and peak stress estimation.

[Figure omitted. See PDF]

Plots of the peak stress data for each confining pressure are used to estimate the bulk strength parameters according to the widely applied Mohr–Coulomb and Hoek–Brown failure criteria (Hoek, 2007; Hoek et al., 2002; Hoek and Brown, 1997) (Fig. 7). The Mohr–Coulomb failure envelopes for the compression tests are shown in Fig. 7a. If tension test results are included, a highly non-linear behaviour of the material is recorded, so that a Mohr–Coulomb envelope is partly not appropriate anymore. Consequently, a non-linear Hoek–Brown envelope is included in Fig. 7b, although it only fits well to all the data for the lacustrine mud material and to low-confining pressure data for the other materials.

The slopes of the elastic parts of the stress–strain curves are used to estimate the bulk elasticity parameters. Figure 7c shows that Young's modulus, $E$, increases with confining pressure (i.e. depth), while Poisson's ratio, $\mathit{\nu }$, shows no trend. Tensile tests reveal lower elastic moduli and Poisson ratios than in the compression tests. An overview of the bulk material properties resulting from these calibration tests is given in Table 4. At low confining pressures, the failure envelopes for bonded particle models are non-linear – cf. Schöpfer et al. (2013). Further details and examples are found in Appendix B4.

Bulk failure envelopes and elasticity parameters of the simulated Dead Sea materials derived from simulated laboratory tests. (a) Mohr–Coulomb failure envelopes, (b) Mohr–Coulomb and Hoek–Brown principal stress, (c) variation of elastic parameters with confining pressure. Left column: low-strength cohesive lacustrine mud; central column: middle-strength alluvial sandy gravel material; right column: higher-strength Holocene salt. CT indicates compression test data; TT indicates tension test data.

[Figure omitted. See PDF]

Bulk material properties of the three investigated Dead Sea materials as derived by simulated rock tests and measurement circles. All values refer to unconfined conditions (i.e. at or close to the surface). Mohr–Coulomb and Hoek–Brown results are based on compression and tension tests on 10 different particle assemblies for each material.

Parameter	Symbol	Unit	Wet lacustrine mud	Alluvial sediment	Holocene salt
Particle packing porosity	$n_{\mathrm{eff}}$	–	0.21	0.2	0.17
Bulk density	${\mathit{\rho }}_{\mathrm{bulk}}$	(kg m${}^{-\mathrm{3}}$)	2145	2200	2075
Young's modulus	$E_{\mathrm{eff}}$	(GPa)	$\mathrm{0.084}\pm {\mathrm{1.2}}^{-\mathrm{2}}$	$\mathrm{0.174}\pm {\mathrm{2.5}}^{-\mathrm{2}}$	$\mathrm{1.106}\pm {\mathrm{126}}^{-\mathrm{3}}$
Poisson's ratio	${\mathit{\nu }}_{\mathrm{eff}}$		$\mathrm{0.19}\pm \mathrm{0.12}$	$\mathrm{0.31}\pm {\mathrm{6}}^{-\mathrm{2}}$	$\mathrm{0.30}\pm \mathrm{0.03}$
Mohr–Coulomb: unconfined compressive	UCS	(MPa)	$-\mathrm{0.25}\pm {\mathrm{5}}^{-\mathrm{3}}$	$-\mathrm{0.52}\pm {\mathrm{8}}^{-\mathrm{3}}$	$-\mathrm{1.23}\pm {\mathrm{1.4}}^{-\mathrm{2}}$
strength
Mohr–Coulomb: unconfined tensile	$T$	(MPa)	$\mathrm{0.2}\pm {\mathrm{4}}^{-\mathrm{3}}$	$\mathrm{0.24}\pm {\mathrm{4}}^{-\mathrm{3}}$	$\mathrm{0.43}\pm {\mathrm{5}}^{-\mathrm{5}}$
strength
Mohr–Coulomb: cohesion	$c$	(MPa)	$\mathrm{0.11}\pm {\mathrm{1.2}}^{-\mathrm{4}}$	$\mathrm{0.18}\pm {\mathrm{5}}^{-\mathrm{4}}$	$\mathrm{0.36}\pm {\mathrm{1}}^{-\mathrm{3}}$
Mohr–Coulomb: friction angle	$\mathit{\varphi }$	(${}^{\circ }$)	$\mathrm{5.7}\pm \mathrm{0.06}$	$\mathrm{22.3}\pm \mathrm{0.17}$	$\mathrm{28.8}\pm \mathrm{0.18}$
Hoek–Brown: unconfined compressive	UCS	(MPa)	$-\mathrm{0.06}\pm \mathrm{0.04}$	$-\mathrm{0.92}\pm \mathrm{0.28}$	$-\mathrm{1.54}\pm \mathrm{0.41}$
strength
Hoek–Brown: unconfined tensile	$T$	(MPa)	$\mathrm{0.01}\pm \mathrm{0.01}$	$\mathrm{0.18}\pm {\mathrm{9}}^{-\mathrm{2}}$	$\mathrm{0.31}\pm \mathrm{0.14}$
strength
Hoek–Brown: ratio compressive/tensile	UCS $/T$	(MPa)	6.0	5.1	5
strength

This calibration shows that the model materials mimic the mechanical responses of the natural materials, and it builds the basis for the analysis of the specific sinkhole formation problem at the Dead Sea, as presented in the following section.

Results of DEM models of void growth and sinkhole collapse

Development in “end-member” Dead Sea materials

We simulated the effect of continuous material removal from a semi-elliptical subrosion zone at 20, 30 or 40 m depth below the initial surface for all three end-member Dead Sea materials. For brevity, we here report on the evolution of the models with subrosion at 30 m depth only; for the detailed evolution of all simulated configurations, see the electronic Appendix.

As shown in Figs. 8 and 9, the evolution of cavity development strongly depends on the mechanical interaction with the surrounding material. The mud is geomechanically the weakest end-member, and even the initial small cavity is not supported by it; the cavity collapses almost instantly after it is generated. Consequently, a cavity of large size (metre scale) never develops in the mud. As material is progressively removed from the subrosion zone, material from around and above the removal zone subsides gradually toward it. A column of subsiding material develops that is partly fault-bound and is characterised internally by downsagging of the overburden layering. This column grows upward until intersecting the surface, where a sag-like sinkhole forms. With further subrosion, the sinkhole grows deeper and wider as areas marginal to the subsiding column slump inwards.

In contrast, the alluvium is strong enough to sustain the cavity as it grows. The growing cavity interacts mechanically with the surrounding material, as sections of the cavity roof and walls collapse into it. Eventually, the overburden above the cavity fails abruptly, and the cavity is closed by the collapse of the overburden into it. The overburden collapse is also usually partly fault-bound with downsagging or with a more complex internal structure. The resultant model sinkhole margins are characterised initially by large and deep (metre-scale) opening mode fractures (ground cracks), inward-tilted blocks and in part by overhanging sides. With further subsidence, the inward-tilted blocks and overhanging sides tend to slump into the sinkhole's centre. The salt is the strongest end-member geomechanically, and so large stable cavities can develop within it – essentially unaffected by deformation of the surrounding material – until only a thin “bridge” of overburden is left.

The mechanical differences in the structural development are highlighted in Fig. 9. For the low-strength mud, stress arching, which tends to stabilise the overburden, is weakly developed around the material removal zone and within the overburden. Stress arching is well developed around and above the cavity in the alluvium, although the absolute values of shear stress are high on the cavity's lateral walls, suggesting that these areas are most susceptible to failure. The stress arch is disrupted upon final failure of the overburden and formation of a sinkhole. In the strong salt, stress arching is best developed and persists even after the thin “bridge” of the remaining overburden fails.

Evolution of DEM model cavity growth and sinkhole collapse in end-member Dead Sea materials. Shown here are selected stages in the development of cavity/sinkhole in salty mud (a), alluvium (b) and rock salt (c). The top row shows the initial cavity growth stage for each material. The layering in the other rows defines passive markers and does not represent any change in material properties.

[Figure omitted. See PDF]

Evolution of maximum shear stress during cavity growth and sinkhole collapse in end-member Dead Sea materials. Shown here are selected stages in the development of cavity/sinkhole in salty mud (a), alluvium (b) and rock salt (c). The same model setups as in Fig. 8 are used.

[Figure omitted. See PDF]

Development in layered Dead Sea materials

We also simulated the effect of continuous material removal from a semi-elliptical subrosion zone at 20, 30 or 40 m depth below the initial surface for layered combinations of the end-member Dead Sea materials. The models comprise a layer of either alluvial sandy gravel or rock salt (0–13 m depth) overlying a lacustrine mud layer (13–40 m depth), followed by the alluvium/salt as a basement, respectively. For brevity, we again report on the evolution of the models with subrosion at a depth of 30 m only (Figs. 10 and 11); for the detailed evolution of all simulated configurations, see the Supplement.

Evolution of DEM model cavity growth and sinkhole collapse in layered configurations of the end-member Dead Sea materials. Shown here are selected stages in the development of cavity/sinkhole in salty mud overlain by alluvium (a), salty mud overlain by rock salt (b) and rock salt overlain by salty mud (c). The top row shows the initial cavity growth stage for each material. Note that the initial cavity again closes rapidly in the mud, leading to a broader zone of subsurface instability.

[Figure omitted. See PDF]

Evolution of maximum shear stress during cavity growth and sinkhole collapse in layered configurations of the end-member Dead Sea materials. Shown here are selected stages in the development of cavity/sinkhole in salty mud overlain by alluvium (a), salty mud overlain by rock salt (b) and rock salt overlain by salty mud (c). The same model setups as in Fig. 10 are used.

[Figure omitted. See PDF]

In general, for layered materials with mud as the subrosion-affected interlayer, the ground tends to fail clearly earlier than for the uniform materials. The mud cannot sustain large cavities and hence fails immediately upon material removal, and the upper mud layers bend. This leads consequently to the development of a cone-shaped underground collapse zone. In alluvium on lacustrine mud, a small subsidence may be noted before collapse and cracks appear even at a certain distance from the main area. Note also the development of ephemeral cavities at the interface with the mud and/or within the alluvium or salt top layers, as deformation migrates upward toward the surface. After the collapse, large and small rotated blocks slump towards the centre, and opening cracks grow downwards to a depth of 12 m around the collapse zone. These blocks define the base of the formed sinkholes. Although salt has double the strength of alluvium, the shapes of the sinkhole for these multilayer models do not differ much, but a small tendency to more overhanging sides is observed. For the condition of lacustrine mud on rock salt, the salt layer sustains large cavity formation, but as soon as the void space reaches $\sim \mathrm{2}$ m below the mud border, the material collapses. The formed sinkhole is a mixture of typical end-member types mentioned in Sect. 1.

As shown in Fig. 11, the mechanical effect of the weak mud layer is to inhibit the development of stable stress arching in the overburden. Where the weaker layer lies below the stronger layer, the development of a collapse zone is indicated as a zone of low stresses, around and above which a stress arching is weakly developed. Effectively, this subsurface collapse zone is mechanically similar to a large cavity. The lack of support from the weak layer concentrates stress in the stronger layer (note the high magnitude of shear stress there), pushing the strong layer toward failure. Where the weak layer overlies the strong layer, the stress arch is well developed until the cavity growth nears the weaker layer. The weaker layer cannot sustain the stress arch, and so the overburden collapses.

Effect of subrosion zone depth

As shown in Fig. 12, the variation of depth of the subrosion zone changes the morphology of the sinkholes. For more details, refer to the electronic Appendix. The removed material in the subrosion zone is assigned a removed “volume” $\mathrm{\Delta }V$, which is based on the area of the removed disk-shaped particles and its unit thickness; see Appendix A.

Sinkhole end-members in dependency of the depth of the subrosion zone for all material combinations investigated in this study. The maximum shear strain is used to visualise the collapse zone.

[Figure omitted. See PDF]

In lacustrine salty mud, for all subrosion depths, the sinkhole collapse is gradual with continuous subsidence. The deeper the subrosion zone, the lower the vertical displacement at the surface, and a greater amount of material needs to be removed before an effect is visible at the surface ($\mathrm{\Delta }V\sim \mathrm{80}$ m${}^{\mathrm{3}}$ for deep vs. $\mathrm{\Delta }V\sim \mathrm{50}$ m${}^{\mathrm{3}}$ for shallow). For a shallow subrosion zone, the sinkholes are $V$ shaped with partly steep margins. For middle subrosion zones, the sinkholes exhibit a compressed $V$ shape with both flat and steep margins. In contrast, the deep subrosion zone leads to bowl-shaped sinkholes with flat sides.

In the homogeneous alluvium models, the sinkhole collapse process varies between sudden (shallow material removal zone) and partly gradual (deep zone). For a shallow subrosion zone, the collapse occurs relatively late at a removed material volume of $\mathrm{\Delta }V\sim \mathrm{400}$ m${}^{\mathrm{3}}$. A long-term stable cavity, also asymmetric (see Fig. 8), can reach the immediate subsurface, and no precursory cracks at the surface appear. The final sinkhole is $A$ shaped with overhanging sides. A deep subrosion zone causes cracking in the overlying layers and at the surface together with subsidence before gradual collapse occurs, commencing relatively early at $\mathrm{\Delta }V\sim \mathrm{80}$ m${}^{\mathrm{3}}$. The final sinkhole is $V$ shaped.

In homogeneous rock salt models, for all subrosion depths, the sinkhole collapse is sudden and occurs after large amounts of material are removed. The sinkholes that form are in all cases $A$ shaped. No surface subsidence can be observed before the collapse, as the void spaces stay stable up to the immediate subsurface. For a shallow subrosion zone, the cavity fails very late at $\mathrm{\Delta }V\sim \mathrm{400}$ m${}^{\mathrm{3}}$, for a middle subrosion zone at $\mathrm{\Delta }V\sim \mathrm{900}$ m${}^{\mathrm{3}}$ and for a deep subrosion zone at $\mathrm{\Delta }V\sim \mathrm{1500}$ m${}^{\mathrm{3}}$. The latter shows pronounced spalling at the sides of the cavity. The shallow model only fails because the material left is of minute thickness.

For the multilayer model alluvium on mud with alluvial basement, the collapse in all cases happens earlier than in pure alluvial material and is sudden. For a shallow subrosion zone, the sinkhole forms at $\mathrm{\Delta }V\sim \mathrm{240}$ m${}^{\mathrm{3}}$; for middle and deep subrosion zones, it forms at $\mathrm{\Delta }V\sim \mathrm{80}$–100 m${}^{\mathrm{3}}$, with little subsidence before collapse onset. For middle and deep subrosion zones, the formed sinkhole is initially narrower but widens with continued material removal. For the shallow zone, this does not happen due to lack of material.

Similar features are observed for the multilayer model salt on mud; the collapse in all cases happens earlier than in pure salt material and is sudden. The removed volume before collapse is similar to results from alluvium on mud; namely, for a shallow subrosion zone the sinkhole forms at $\mathrm{\Delta }V\sim \mathrm{240}$ m${}^{\mathrm{3}}$; for middle and deep subrosion zones, it forms at $\mathrm{\Delta }V\sim \mathrm{80}$–120 m${}^{\mathrm{3}}$, with little pre-collapse subsidence and compression ridges. The sinkhole morphologies are similar to the ones for alluvium on mud, but a tendency to larger block size and a more pronounced overhanging is observed.

Thickness-to-diameter ratio at the onset of collapse

Figure 13a shows the estimated $T/D$ ratios at the onset of the collapse ($T/D_{\mathrm{crit}}$) for all model setups independent of the subrosion zone depth. A collapse hereby is defined when both particle movement at the surface and in the subsurface occurs. A distinction for the different involved materials is found. Pure lacustrine mud models generally fail at higher ratios of $T/D_{\mathrm{crit}}\ge \mathrm{0.5}$ than the majority of the other models. Multilayer models with mud underlying alluvium or salt show low $T/D_{\mathrm{crit}}\le \mathrm{0.5}$, while pure alluvium or salt models have the highest and lowest measured values, respectively. A collection of the mean values is given in Table 5. The deeper the subrosion zone in both multilayer and uniform material models, the less material needs to be removed to trigger a collapse (Fig. 13b).

Critical thickness to diameter ratios for modelled sinkhole collapse onsets. The error is based on the mean between different particle assemblies for each setting.

Subrosion zone depth/	Lacustrine mud	Alluvium	Salt	Alluvium on mud	Salt on mud	Mud on salt
modelled material setup
Shallow	$\mathrm{0.66}\pm \mathrm{0.01}$	$\mathrm{0.06}\pm \mathrm{0.05}$	$\mathrm{0.02}\pm \mathrm{0.01}$	$\mathrm{0.22}\pm \mathrm{0.05}$	$\mathrm{0.17}\pm \mathrm{0.1}$	–
Middle	$\mathrm{0.57}\pm \mathrm{0.05}$	$\mathrm{0.5}\pm \mathrm{0.29}$	$\mathrm{0.03}\pm \mathrm{0.03}$	$\mathrm{0.32}\pm \mathrm{0.08}$	$\mathrm{0.48}\pm \mathrm{0.15}$	$\mathrm{0.43}\pm \mathrm{0.1}$
Deep	$\mathrm{0.57}\pm \mathrm{0.02}$	$\mathrm{0.97}\pm \mathrm{0.1}$	$\mathrm{0.08}\pm \mathrm{0.01}$	$\mathrm{0.27}\pm \mathrm{0.06}$	$\mathrm{0.51}\pm \mathrm{0.03}$	–

Parameters at the onset of collapse. (a) Thickness to diameter ratio vs. removed volume at the onset of collapse. Note that the $T/D_{\mathrm{crit}}$ ratio is either referring to the (stable) subsurface cavity or the (instable) subsurface collapse zone. Red colours stand for mud, orange for alluvial material and blue for salt and the respective multilayer models. (b) Removed volume at the onset of collapse in dependency of the subrosion zone depth.

[Figure omitted. See PDF]

Comparison with data derived from photogrammetry at Ghor Al-Haditha sinkhole site

Surface displacement

In Fig. 14, we compare the topographic profiles of sinkholes derived from photogrammetric studies at Ghor Al-Haditha (see Sect. 1, Fig. 3) with our results from DEM sinkhole modelling. In Fig. 14a, we show the simulated sinkhole morphologies for different evolution stages for a subrosion zone with intermediate depth (30 m). To facilitate the comparison, the topographic profiles derived by photogrammetry have been normalised and the axes have been adjusted to the same dimensions as for the models (Fig. 14b). An impressive similarity can be found for these sinkhole end-members both in terms of lateral extent and subsidence amplitude: (1) the mud sinkhole in the field appears to be of an early-stage sinkhole but with a larger extension laterally; (2) the alluvium sinkhole shape is remarkably similar to the late-stage (evolved) modelled sinkholes both laterally and vertically; (3) the salt sinkhole is comparable to the respective simulation result for an early-stage salt sinkhole.

Topographic cross-sections of sinkholes in different cover materials. (a) Modelled profiles for subrosion zones in 30 m depth and three different evolution stages. Plotted is the vertical position of the surface particles for five different particle assemblies of each tested model setup. Model sets are pure lacustrine mud (left row), alluvium on mud multilayer (middle row) and salt on mud multilayer (right row). (b) Field data replotted topographic profiles of the three different sinkhole morphologies from Fig. 3. The distance and altitude are normalised for better comparison with the models.

[Figure omitted. See PDF]

These findings are essentially confirmed by knowledge about the rather recent development of the sinkholes selected in the mud and salt flats and the older, more evolved sinkholes in the alluvial fan of Ghor Al-Haditha (Al-Halbouni et al., 2017). Our models, which are based on realistic material parameter estimation, hence reproduce the topographic features of the sinkholes successfully in the field site. This result is even better reflected in the De $/$ Di analysis described in the following section.

Sinkhole depth $/$ diameter ratios

In Fig. 15, we compare the sinkhole depth $/$ diameter (De $/$ Di) ratios for the DEM models and for natural equivalents at the Ghor Al-Haditha field site. We use results from the models with uniform lacustrine mud from the layered models of alluvium on mud for different stages of collapse (early, middle and late). Figure 15 does not include data for salt-on-mud layered models, as natural equivalents were not mapped by Al-Halbouni et al. (2017), but the De $/$ Di ratios for these models are given in Table 6 and are generally similar to those in alluvium-on-mud models. Mean values of De $/$ Di are $\mathrm{0.37}\pm \mathrm{0.15}$ for alluvium, $\mathrm{0.15}\pm \mathrm{0.02}$ for mud and $\mathrm{0.33}\pm \mathrm{0.11}$ for salt, and close to the statistical estimates given by Al-Halbouni et al. (2017) and to the examples shown in Fig. 2.

In simulations with uniform mud material, the fit to the De $/$ Di data is good in the early and intermediate stages of collapse (Fig. 15a). For the late stages of these models, the model De $/$ Di ratios are at the outer bound of the natural range. Conversely, in simulations with alluvium-on-mud material, the fit to the De $/$ Di data is better in the intermediate and early stages of collapse (Fig. 15b). For the early stages in alluvium, the results are at the lower margin of the depth.

Sinkhole depth $/$ diameter ratios from photogrammetry (Al-Halbouni et al., 2017) and DEM model results of this study. For field data, depth and diameter for the materials alluvium and mud were determined for 237 sinkholes. The mean values are $\mathrm{0.4}\pm \mathrm{0.11}$ for sinkholes in alluvium and $\mathrm{0.14}\pm \mathrm{0.04}$ for those in mud. No field data were available for salt at the time of this study. For the models, we distinguish between early collapses (circles), middle collapses (triangles) and late collapses (squares); (a) mud-flat sediments and lacustrine mud models, (b) alluvial sediments and multilayer alluvium on lacustrine mud model results.

[Figure omitted. See PDF]

Depth $/$ diameter ratios for modelled sinkhole collapses. A set of five models for each material combination is analysed for a subrosion zone at middle depth (30 m). De $/$ Di ratios for alluvium and salt are generally higher than for mud.

Collapse stage	Removed volume (m${}^{\mathrm{3}}$)	Lacustrine mud	Alluvium on mud	Salt on mud
Early	100	$\mathrm{0.09}\pm \mathrm{0.01}$	$\mathrm{0.13}\pm \mathrm{0.06}$	$\mathrm{0.08}\pm \mathrm{0.05}$
Middle	160	$\mathrm{0.13}\pm \mathrm{0.02}$	$\mathrm{0.25}\pm \mathrm{0.1}$	$\mathrm{0.21}\pm \mathrm{0.08}$
Late	220	$\mathrm{0.15}\pm \mathrm{0.02}$	$\mathrm{0.37}\pm \mathrm{0.15}$	$\mathrm{0.33}\pm \mathrm{0.11}$

Discussion

Comparison to previous DEM and non-DEM studies of cavity generation and sinkhole collapse

Baryakh et al. (2008, 2009) used the DEM to investigate the effect of depth, geometry and mechanical properties on the collapsed state in karst. Their approach is to some extent similar to ours; however, essentially only the position of a rectangular or an arched cavity was varied for different uncalibrated materials. In contrast, our numerical simulations allow for a mechanical interaction of a slowly growing void space with the surrounding rock and provide calibrated bulk rock parameters. Consequently, the material removal either creates a cavity or not, leading to variably shaped subsurface collapse zones, details of which are elaborated on later. Hatzor et al. (2010) used jointed blocky rock mass (DFN) modelling to define stability criteria ($T/D$ ratios) for a rectangular cavern in high-strength (UCS > 10 MPa) rocks. One conclusion of their study, namely the conservative $T/D$ $=$ 1.0 for large cavity sizes, may also apply to our results for homogeneous, relatively weak rock and cohesive soil models. Nevertheless, the stability depends strongly on the collapse zone geometry, and the well-known stability limit for deep-seated excavation from Terzaghi (1946) does not hold for our shallow collapse zones. A FEM approach from Shalev and Lyakhovsky (2012) addresses sinkhole formation by utilising viscoelastic rheology with a damage model. It is applied to the sinkhole hazard at the Dead Sea and relates the different deformation modes (viscoelastic vs. brittle) to the different mechanical properties of the involved materials (mud vs. alluvium) and their common morphological characteristics. However, no field data comparison is given, and sinkhole formation is only simulated using a simplified cavity geometry that does not evolve.

In summary, earlier studies lack a detailed calibration of the model strength parameters to field and laboratory estimates, and quantitative comparisons of model results with measured data are limited or absent. Our study hence fills this important gap and explicitly simulates cavity growth and related sinkhole development and therefore provides a significant advance in this field.

Model testing, and benchmarking and limitations

Our tests and model benchmarking provide several new insights for undertaking the simulation of karstic void development and sinkhole collapse under gravity with the DEM. As expected, there is a strong sensitivity of model results (displacement) not only to parameters such as model dimensions and resolution but also to model shape, with the best results attained for relatively high-resolution and equidimensional model setups. Our tests also show that the method of cavity generation has only a minor impact on the surface displacement pattern. Cavity generation by particle deletion differs from generation by particle radius reduction mainly in the much longer model runtime for the latter. This is reasonable given the elastic and quasi-static conditions of the DEM test models. By such tests, we infer that the models with non-elastic deformation (i.e. cavity wall failure and sinkhole collapse) are also insensitive to cavity generation method as long as they are run under quasi-static conditions, as was the case in our study.

In the benchmarking tests, the DEM surface displacements for a circular cavity in a gravitationally loaded elastic material closely resemble those predicted by the BEM model and the Kirsch solution both in the far and near fields of the subsidence centre (see Fig. 5). A perfect match is not expected, despite our efforts to compare like for like, having in mind the intrinsic differences between these models in terms of material properties and boundary conditions. The Kirsch results nonetheless provide the best match to the DEM results for both vertical displacement and displacement differences. Overall, the DEM and Kirsch curves fit in the near field and behave similarly and realistically (tendency to zero) in the far field. The BEM models offer a plane-strain solution for a hydrostatic remote stress, while the two-dimensional DEM model does not consider out-of-plane stress–strain and additionally has before cavity creation a horizontal to vertical stress ratio $K_{\mathrm{0}}$ equal to $\sim \mathrm{1}/\mathrm{4}$. This leads to the generally narrower vertical and horizontal displacement curves in the DEM models at the centre of the subsidence area. The mismatch to the Mindlin solution is greatest in the far-field displacements; these displacements as predicted by the Mindlin solution seem unrealistic given that they progressively increase away from the cavity. Consequently, for the purposes of this work and in light of the minute differences between the DEM results and the BEM/analytical solutions (sub-millimetre for displacements and micrometer for displacement differences, except for the Mindlin solution in the far field), we consider the DEM model approach here to be a valid numerical approximation of the problem.

The manner of cavity growth and its timing relative to collapse are, of course, simplified approximations to complex processes of dissolution and mechanical erosion of the subsurface as they occur in nature. The model cavity grows by instantaneous and repeated material removal of the same volume within a domain of simplified shape. In reality, cavity growth may occur on extremely long to relatively short timescales, depending on the nature of the materials (e.g. limestone vs. salt) and hydrogeological conditions (e.g. porous flow, conduit flow, dripping, flash floods). The cycling to quasi-static equilibrium during each model growth increment ensures, however, that cavity growth rate is smaller than or equal to collapse rate, as expected in nature. An improvement will be to adjust the cavity area growth function to follow typical dissolution laws (cf. Dreybrodt and Kaufmann, 2007; Kaufmann and Romanov, 2016) and thus to develop more complex and realistic cavity geometries.

Geomechanical parameter calibration

The outcomes of the simulated compression and tension tests (Table 4) closely agree with literature values and estimations from geotechnical studies and seismic velocity measurements (Table 2), in terms of UCS ranges, bulk densities, Young's modulus and Poisson ratios. The friction angles of the simulated sand–gravel and rock salt materials are slightly lower than the desired values but fit well in the case of the low-strength lacustrine clay material. Low-friction angles are typical for bonded particle models (cf., e.g. Schöpfer et al., 2017), because both sliding and rotation of particles accommodate bulk deformation; with the contact model used in the present study, the latter cannot be inhibited even with large friction coefficients. It is well known from other DEM studies that UCS $/T$ ratios in bonded-particle materials are lower when compared to natural rock (UCS $/T\sim \mathrm{10}$) and soils (UCS $/T\sim \mathrm{8}$; Koolen and Vaandrager, 1984), reflecting the discretisation by means of circular/spherical particles (Schöpfer et al., 2007, 2009).

The Mohr–Coulomb and Hoek–Brown failure envelopes (Fig. 7) fitted to the calibration data serve as guides to the material behaviour. These envelopes were chosen as they are widely used in geomechanics, and so the overall behaviour of the model materials is readily assessed from them. In detail, however, neither envelope provides a perfect fit to the calibration results. This may be a consequence of the timing of confinement of the particle assembly, which is here done before installing the parallel bonds. Our results indicate that this may lead to some stress-path-dependent behaviour that is more complex than can be represented fully by either Mohr–Coulomb or Hoek–Brown envelopes. A thorough exploration of such complexities is well beyond the scope of this paper, but it could be subject of future work.

It is well known that the relationship between field-scale rock parameters and those determined at the laboratory sample scale depends strongly on the degree of fracturing or alteration of the rock mass (Schultz, 1996). Given that the materials we studied are of rather low strength and are weakly consolidated materials (in contrast to hard karst rock in which sinkholes often form), we neglected the effect of pre-existing weaknesses (e.g. tectonic fractures). We hence adopted literature values for salt and mud derived from laboratory-scale measurements. A poorly understood effect in the Dead Sea materials is, however, the influence of water content which may lead to time-dependent geomechanical behaviours (see Shalev and Lyakhovsky, 2012) that is not accounted for in our models. In principle, however, the modelling scheme we developed could be adapted to account for time-dependent (e.g. viscoelastic) material behaviour.

General implications for cavity and sinkhole formation

Structural or morphological features of sinkholes

The DEM models of sinkhole collapse show a wide range of structural or morphological features that are found at natural sinkholes, and they highlight how these features reflect the mechanical properties of the material in which the sinkholes form. Similar near-surface structural features are found at volcanic collapse calderas and pit craters, and similar explanation in terms of mechanical properties of the near-surface materials have been proposed (Holohan et al., 2011; Poppe et al., 2015).

In relatively weak materials (here the simulated “mud”), the near-surface strain is distributed across many small fractures, such that there is no sharp margin to the sinkhole. Subsidence at the surface develops gradually before the collapse develops (if at all) and the material's response is brittle–ductile. The sinkhole also widens gradually as it deepens. Overall, the sinkhole formation process is similar to classic “cover sagging” or “cover collapse” with partial suffusion (cf. Gutiérrez et al., 2008).

In relatively strong materials (here the simulated “alluvial sand–gravel” and “salt”), the strain is localised on fewer but larger fractures that develop as faults (shear fractures) and/or deep cracks (opening-mode fractures). Structures like compression ridges form in the centre of the subsidence area. Sinkhole margins in such materials are consequently sharp, steep and, at least initially, overhanging. Any subsidence before collapse is slight, although this depends partly on material rigidity (i.e. modulus); the material's response is brittle. The sinkhole also widens as it deepens but in more of a stepwise manner as new marginal fractures form and delimit marginal blocks. Overall, the collapse style is similar to classic “caprock collapse” or “bedrock collapse” (see Gutiérrez et al., 2008). In extremely strong materials, there may be little or no collapse at all – in the limit, the hole may result simply from the intersection of an essentially stable, growing cavity with the ground surface.

Stability of cavities and relationship to sinkhole geometry

The stability of cavities in the DEM models is clearly related to the strength of the material and to the depth of the material removal zone. In general, the cavity stability depends on a combination of material strength (UCS, $T$, friction coefficient) and geometric properties (cavity geometry, $T/D$ ratio). In principle, larger $T/D$ and stronger materials promote larger void spaces in the underground as stable compression arches build up (Fig. 16; see also Holohan et al., 2015). Thus, for a given $T/D$, cavities are unstable in the weak “mud” material but are stable in the stronger “sand–gravel” and “rock salt” materials. As the cavity grows, however, the $T/D$ ratio decreases and ultimately the overburden geometry can no longer support its weight. Eventually, the overburden will fail partially or completely and collapse into the cavity.

Maximum compressive stress for representative models after the same amount of material removal ($\sim \mathrm{33}$ m${}^{\mathrm{3}}$). Different depths of subrosion zones are compared: (a) shallow (20 m), (b) middle (30 m) and (c) deep (40 m). The deeper the zone, the higher the maximum compressive stress above the created void space. The stronger the material, the more pronounced the compression arch. Tensile stresses are observed directly above the removal zone.

[Figure omitted. See PDF]

The gravitational stress field in the models also means that the absolute depth, and not just relative depth as expressed by $T/D$, is critical, however. The deeper the cavity, inside which stresses are zero, the higher the differential stress immediately around it (Figs. 9, 11 and 16). This accounts for the observation in our models that, for a given material strength, deeper-seated cavities fail earlier than shallow ones in these weakly consolidated materials. Overall, our results indicate that cavity sizes and stability, and hence the style of sinkhole collapse, will depend on the material strength and depth of dissolution. Thus, caprock collapse sinkholes, which form above large cavities (Fig. 1), may be favoured for relatively strong material and/or shallow dissolution levels. Dropout or suffusion sinkholes may be favoured by relatively weak material and/or deep dissolution levels. In the limit, no macro-scale cavities will form below a certain dissolution zone depth in a given material, as in situ stresses become too high for that material to support such cavities.

The DEM models also show how the interaction of material removal and mechanical instability can lead to cavity growth. This is seen mainly in moderately strong DEM material (here the “sand and gravel”), where void spaces usually stay stable until large volumes of material are removed, with typical spalling at the sides rather than from the roof (Fig. 8). This lateral spalling of the cavity is typical of “tunnel breakouts” encountered by engineers and arises from the in situ stress field in the DEM model surrounding the cavity being characterised by a $K_{\mathrm{0}}$ < 1 (${\mathit{\sigma }}_{xx}$ < ${\mathit{\sigma }}_{yy}$). In nature, a feedback mechanism may arise from such spalling, whereby lateral or vertical spalling expose more fresh surfaces to dissolution and reduces the overburden $T/D$, leading to further cavity growth and instability, leading to more spalling, etc.

Another important result of our DEM models is that multilayer models with a weak (mud) interlayer fail earlier than the models with a uniform material. This is not only because the integrated strength of the overburden is lessened, but also because the rapid failure of any cavities in the weak layer effectively increases the stress concentration in the strong overlying layer, similar to a beam (Fig. 16), leading to bending induced stresses with inner arc contraction and outer arc extension. This is contrary to the higher $T/D$ ratios for the same amount of removed volume in the homogeneous layer models in which a stable cavity develops.

A consequence of such material-controlled cavity stability is that, as is often inferred for nature (e.g. Waltham et al., 2005), the geometric relationship between subsurface cavities and sinkholes is not a straightforward one. In the weak DEM model material, a sinkhole can have little or no geometric relationship to a cavity, because cavities are not sustained at any comparable scale. In the strong DEM model materials, on the other hand, the sinkhole geometry may relate to cavity geometry to a variable degree. This relationship may be especially direct in the case of a shallow removal zone and a very strong material, where a cavity can stably grow upward with little or no collapse until intersecting the ground surface. Overall, our results reinforce the point that the use of continuum-based methods to estimate cavity geometry from sinkhole geometry (i.e. where there are large permanent strains) should be treated with caution (see also Fuenkajorn and Archeeploha, 2010 and Holohan et al., 2017).

Future work will include a variation of lateral (long-wall-mining-like), vertical (tube-like) and multiple void space growth systems. Especially for typical karst simulations, multiple void spaces with different growth functions and geometries are a more suitable, complex approach. Another aspect is the role of hydrostatic (buoyancy) and pore pressure, which is usually an important factor regarding soil liquefaction and landslides due to the reduction of effective stress (cf., e.g. Tharp, 1999; Zeev et al., 2017; Clément et al., 2018) and has been ignored in these simulations for simplicity. A possible DEM approach is to apply forces to the boundary particles of the void space to simulate the pressure inside a water-filled cavity or to apply forces related to the pore spaces between particles to simulate hydrofractures (Yoon et al., 2015). An alternative is the combination of FEM and DEM with accounting for drag forces due to fluid flow or other combined particle-lattice model schemes (Ghani et al., 2013).

Implications for sinkhole formation at the Dead Sea

In general, the good fit of model sinkhole geometry with the observed topography of sinkholes at Ghor Al-Haditha (Sect. 4) confirms the suitability of the DEM approach and allows for interpretation of morphological features there. In addition, structures as found in the simulations are visible also in the field, such as sagging layers and distributed marginal fracturing in weak materials, as well as cavities, compression ridges (pop-up structures) and overhanging sides in stronger materials. For a still better fit to the low diameter results of the field (Fig. 15), we would need to use a wider variety of the void space growth functions, geometries and subrosion zone depths, as expected to happen in nature. Due to computational costs, this has not been included in this study. Nonetheless, the already good agreement between the paths of depth $/$ diameter of the existing model sinkholes as subsidence evolves and the distribution of depth $/$ diameter values in the field (Fig. 15) strongly suggests that those distributions represent growth trends of the natural sinkholes that are controlled ultimately by material properties (Fig. 17).

Simulated sinkhole depth $/$ diameter interpretation. The simulations reveal a tendency towards deeper sinkholes in alluvium and both deeper and wider sinkholes in mud, a trend that is able to explain the observations for sinkholes at Ghor Al-Haditha.

[Figure omitted. See PDF]

Since material heterogeneity is the rule rather than the exception in nature, and since our simulation results fit well to seismic and photogrammetric studies in the area of Ghor Al-Haditha (Al-Halbouni et al., 2017; Polom et al., 2018), we consider our multilayer models as favourable over uniform models for Ghor Al-Haditha. The exact values of large-scale material strength, however, due to the described material testing procedure with a constant particle packing porosity and the limitations of literature laboratory scale values under the assumption of intact rock, should be rather used carefully. Lower strength for the materials in the field is highly probable, as the observed maximum crack depth in alluvial and salt materials (4 m) is less than in the DEM simulations (up to 12 m). This is probably because in detail “pure” sand–gravel or rock salt is rare on a large scale at the site – usually, there is plenty of muddy material interbedded. However, some general observations for the models with materials and material successions typical at the field site of Ghor Al-Haditha in Jordan can be drawn from the following simulations:

• A weak lacustrine mud layer beneath a strong cover material favours sinkhole formation. Even high-strength material like the salt would collapse in such a setting.

• A middle–deep subrosion zone (30–40 m) leads to collapses even for the pure alluvium models, which means that a subrosion acting only in the alluvial sediments can similarly cause sinkhole formations like those with a weak interlayer. Only a higher-volume removal is needed.

• The pure salt models do not produce typical sinkholes as observed in the field zone. This fact can be related either to a lack of such a thick and strong cover material in nature or a too-high strength assigned in the model. It is perhaps worth noting that at the Lisan Peninsula, close to the field area at the Dead Sea, large (several-metre-scale) cavities and arches were observed here in Holocene Dead Sea salt (Closson et al., 2007). On the other hand, the observed salt exposure at our field site contains rather thin salt layers, interleaved with mud on a centimetre scale, so that the bulk material strength there is expected to be lower than that simulated.

• The possibility to record surface subsidence before actual sinkhole collapse depends on both the cover material type and the depth of the subrosion zone. A multilayer model of a middle-deep subrosion zone with a large subsurface collapse zone may produce recordable surface signatures in the order of sub cm before the onset of collapse.

Finally, the single void collapse concept explored in this paper may sufficiently explain some individual sinkhole occurrences at Ghor Al-Haditha and elsewhere around the Dead Sea (cf. laboratory experiments by Oz et al., 2016); the coalescence, sequence evolution and sinkhole cluster structures, morphological expressions at the surface and larger sinkhole depression areas may not. For this, a more sophisticated approach of multiple void space growth, testing different geometries and a more realistic subrosion process is necessary and will be addressed in a future paper.

Summary and conclusions

In this work, we presented a benchmarked and calibrated 2-D distinct element modelling approach to simulating the process of both cavity growth and sinkhole development. Our principal findings are as follows.

Firstly, we presented a computationally fast DEM approach to simulating sinkhole formation by instantaneous, quasi-static, stepwise material removal in a single void space at a depth of an arbitrarily shaped geometry under gravitational loading. We successfully benchmarked the models with analytical and BEM solutions yielding a sub-millimetre degree of agreement for surface displacements and displacement differences.

Secondly, we performed simulated compression and tension tests to determine microscopic bond strength parameters and moduli calibrated by intact rock literature values and field estimates for the three materials common at the Dead Sea shoreline. The simulated rock tests yield low bulk strength (UCS $\sim \mathrm{0.06}$–0.25 MPa) for lacustrine mud, middle bulk strength (UCS $\sim \mathrm{0.53}$–0.92 MPa) for alluvial sandy gravel sediments and high bulk strength (UCS $\sim \mathrm{1.23}$–1.54 MPa) for rock salt materials, based on Mohr–Coulomb and Hoek–Brown fits.

Thirdly, we simulated a cavity growth until sinkhole collapse in uniform materials. Cavity development is controlled by the interaction of the material strength and the depth of material removal. Weak materials do not support large cavities, and so subsidence is characterised by gradual sagging and suffusion-type collapse into the material removal zone. Stronger materials support the development of large cavities at the material removal zone, leading to sinkhole formation by the sudden collapse of the overburden (caprock or cover collapse type sinkholes). At one end of the spectrum, near the Earth's surface, very strong materials may support cavity growth until the intersection with the ground surface, giving rise to sinkholes with little or no collapse. At the other end of the spectrum, below sufficient depth and for a given material strength, the development of cavities on a significant scale is inhibited as gravitational stresses are too high.

Fourthly, we simulated a cavity growth until the sinkhole collapse in multilayered materials. We show with the inclusion of weak layers, either as cover material or as subroded bedrock material, results in sinkhole development with less volume of removed material than in the case of uniform model material. Such development is not only due to an integrated weakening of the overburden but also due to the growth of a subsurface collapse zone in the weak material that geometrically destabilises the overburden.

Lastly, we compare the developed morphologies from a set of models for all three materials with photogrammetric analysis from the sinkhole area of Ghor Al-Haditha in Jordan. Our approach produces physically realistic sinkhole shapes and successfully reproduces typically measured sinkhole depth $/$ diameter ratios of 0.15 in mud-flat material, 0.37 for sinkholes in alluvium and 0.33 in salt. The field distribution appears hereby to be related to evolution stages of the sinkholes between early and late collapses. A weak (mud) interlayer and/or a deeper lying subrosion zone enhances formation of sinkholes in materials typical of the Dead Sea margins.

The distinct element method and its use for simulating geomaterials

The DEM is a specific scheme of undeformable particles and deformable contacts developed by (Cundall, 1971). In the PFC2D v5.035 software, the DEM is used to implement Newton–Euler equations of motion and rotation on disk-shaped particles (Itasca Cooperation Group, 2014; Potyondy, 2014a; Potyondy and Cundall, 2004) (Fig. A1a). In PFC, the resolution scheme is an explicit second-order velocity Verlet algorithm (Verlet, 1967). The particles are assigned a mass and a radius, are initially unbonded and are free to move and rotate depending on external forces. Particles interact only at contact points between particles and wall facets, where the mechanical interaction is treated in terms of a frictional contact with a set of linear elastic springs that are assigned normal and shear stiffness (Fig. A1b). The “rigidity” of the particles is defined by setting the elastic Young constant in accordance to the spring stiffness. An additional bonding of the elements can be performed, whereby many different bond types can be specified. Here, we use the parallel-bond model (Potyondy and Cundall, 2004), which is defined in terms of a set of linear elastic springs in parallel to the linear contact bond. The parallel bonds allow for tensile forces and bending moments between the bonded particles, and they break once their strength is exceeded. Here, we set the bonds to have the same material constants (microproperties) as the particles, like stiffness and elastic modulus, but since bond strength is defined similar to a Mohr–Coulomb failure criterion, the bonds are also assigned a cohesion, tensile strength and friction angle (Fig. A1c).

Schematic description of 2-D DEM modelling with PFC2D v5. (a) Inter-particle and particle-wall force chains developed after gravity settling of an assembly of balls in a box of dimension $H\times W$. The close-ups below show the pore space, contact planes and nomenclature. Particles, although undeformable, are allowed to overlap slightly or have a small gap $g_c$. In both cases, linear contacts and optional parallel bonds are active. These bonds act additionally to the linear contacts. (b) Close-up of the notional and contact planes with all elements necessary for the physical definition of the contact and bond interactions. (c) Failure criterion for parallel bonds. Compression in this study is considered as negative.

[Figure omitted. See PDF]

Mathematical details of the DEM method implemented in the PFC software

The Newton–Euler equations are solved in a finite difference explicit time-stepping algorithm involving dynamic relaxation (Cundall, 1971; Jing and Stephansson, 2007). During the procedure, Newton's second law and the force-displacement law is solved for each of the particles and its contacts (Potyondy and Cundall, 2004). For a 2-D system of coupled rigid elements, the differential equations solved by the explicit time-marching relaxation scheme for a particle of mass $m$ are (Jing and Stephansson, 2007) $$\begin{array}{cc}{\displaystyle }m{\ddot{u}}_x^t+\mathit{\alpha }m{\dot{u}}_x^t& {\displaystyle }=F_{x}m{\ddot{u}}_y^t+\mathit{\alpha }m{\dot{u}}_y^t=F_y,\\ & {\displaystyle }& {\displaystyle }I{\ddot{\mathit{\theta }}}^t+\mathit{\alpha }I{\dot{\mathit{\theta }}}^t=M,\end{array}$$ with $F$ as force, $u$ as displacement, $\dot{u}$ as velocity, $\ddot{u}$ as acceleration, $\mathit{\alpha }$ as damping, $M$ as moment of force, $I$ as inertia, $\mathit{\theta }$ as Euler rotation angle, $\dot{\mathit{\theta }}$ as Euler rotation velocity and $\ddot{\mathit{\theta }}$ as Euler rotation acceleration at a certain time $t$.

It is assumed that (1) velocities and accelerations within one time step are constant and (2) that the step chosen is small enough that disturbances, which occur due to external or body forces, particle or boundary wall movement, propagate only to the neighbours of the particles. The resulting velocity and acceleration components for both the translational and rotational motion of one particle are determined via a finite difference scheme successively for each time step $t$ (Jing and Stephansson, 2007): $$\begin{array}{cc}{\displaystyle }{\ddot{u}}_i^t& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{{\left(\mathrm{\Delta }t\right)}^{\mathrm{2}}}}\left(u_i^{t+\mathrm{1}}-\mathrm{2}{u}_i^t+u_i^{t-\mathrm{1}}\right)\\ & {\displaystyle }& {\displaystyle }{\dot{u}}_i^t={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }t}}\left(u_i^{t+\mathrm{1}}-u_i^{t-\mathrm{1}}\right),\end{array}$$ with $i$ as $x$ or $y$ and the equations for $\dot{\mathit{\theta }}$ and $\ddot{\mathit{\theta }}$ accordingly. The displacement calculation is generally one time step ahead of velocities'/accelerations' calculation, and constitutive laws of arbitrary complexity (Jing and Stephansson, 2007) can be added between the contacts without numerical instability. The kinematic critical time step $\mathrm{\Delta }{t}_{\mathrm{crit}}=\min \left(\sqrt{\frac{m_i}{k_i}}\right)$ is determined for an infinite multiple set of masses $m_i$ and springs with stiffness $k_i$ to allow for the above constraints and solution of the equations.

The equilibrium is defined by a convergence criterion, where the ratio between the “out-of-balance” forces to the overall forces is below a defined threshold (solve ratio, SR) of usually ${\mathrm{1}}^{-\mathrm{5}}$ or lower. This “solving” can be performed for the mean (SR${}_{\mathrm{mean}}$) or maximum (SR${}_{\max }$) forces that appear in the model. A problem can occur when absolute normal force calculation during material gravity settling is used: interlocked forces due to the particle overlap may not be released during further simulation. This issue has been overcome by introducing an incremental normal force calculation (Fakhimi, 2004) which is implemented in PFC2D v5 (Itasca Cooperation Group, 2014; Potyondy, 2014b).

Creating a gravitationally loaded synthetic rock mass in the DEM

Creation of a bonded particle assembly in this study followed that of Holohan et al. (2011) and involved the following chain of steps:

• Creation of an unbonded particle assembly of defined particle sizes, porosity and geometrical distribution. A uniform distribution of particles between a defined minimum and maximum radius is placed randomly in the model box of size $H\times W$. The unbonded material is limited by three walls with low-friction (0.01) elastic interaction. The radii distribution in this study is chosen to be equal between the minimum and maximum assigned radius (Table B1) according to the defined porosity.

• Assignment of material domains. The mechanical properties are distributed in the assembly according to the desired model setup, e.g. layering. The linear contact model is installed between each two entities at a distance smaller than or equal to the surface gap.

• Gravity settling. Gravity acts as the main body force. A settling criterion is applied; i.e. the material is considered as settled when a certain threshold, here SR${}_{\mathrm{mean}}={\mathrm{1}}^{-\mathrm{6}}$, of the velocity and displacement change of the particles between two time steps is reached. The material is settled under low friction until the defined solve ratio.

• Particle bonding. The created assembly contains, as real rock, interlocked forces. Now the particle bonding is applied according to a chosen bond type (parallel) and the model is cycled into equilibrium. Linear contact friction is set to the defined value.

At each step, the material assembly is cycled until a static equilibrium is reached. The behaviour of a DEM model depends strongly on the material packing assembly (Schöpfer et al., 2009), and so a spectrum of solutions is usually obtained by performing multiple realisations for different assemblies. The above chain is thus repeated to produce many random particle assemblies that may be used to obtain a statistical mean of packing-dependent model outcomes. In this study, the procedure was repeated generally for 5–10 random assemblies of the particles.

Details on model convergence tests, benchmarking, material calibration and the final sinkhole model implementation

The following section gives an overview over the performed DEM model convergence, void space installation and benchmarking tests that were performed to determine the optimal sinkhole formation modelling setup. Table B1 summarises the main DEM model parameters used for the tests.

Dimensions of the model and contact and particle properties used in development and testing of DEM cavity formation models.

Geometric parameter or microparameter	Symbol	Unit	Value range (common)
Model height	$H$	(m)	100–500 (400)
Model width	$W$	(m)	100–500 (400)
Material porosity	$n$		0.2
Minimum particle radius	$R_{\min }$	(m)	4.98–0.17 (0.24)
Mean particle radius	$\stackrel{\mathrm{‾}}{R}$	(m)	6.65–0.23 (0.32)
Particle radius factor	$R_{\mathrm{f}}$		1.66
Particle density	$\mathit{\rho }$	(kg m${}^{-\mathrm{3}}$)	2500
Solve ratio unbalanced/balanced forces	SR		${\mathrm{1}}^{-\mathrm{6}}$
Cavity centre depth	$h$	(m)	35
Cavity radius	$r$	(m)	5
Boundary walls Young's modulus	$E_{\mathrm{W}}$	(GPa)	$-\mathrm{5}$
Parallel bond Young's modulus	$E_{\mathrm{Pb}}$	(GPa)	$-\mathrm{5}$
Parallel bond tensile strength	${\stackrel{\mathrm{‾}}{\mathit{\sigma }}}_c$	(MPa)	1000
Parallel bond cohesion	$\stackrel{\mathrm{‾}}{c}$	(MPa)	1000
Parallel bond friction angle	$\mathit{\varphi }$	(${}^{\circ }$)	30
Parallel bond ratio normal/shear stiffness	${\stackrel{\mathrm{‾}}{k}}_{\mathrm{n}}/{\stackrel{\mathrm{‾}}{k}}_{\mathrm{s}}$		2.5
Linear contact Young's modulus	E	(GPa)	$-\mathrm{5}$
Linear contact friction coefficient	$\mathit{\mu }$		0.01–0.5 (0.5)
Linear contact normal/shear stiffness ratio	$k_{\mathrm{n}}/k_{\mathrm{s}}$		2.5
Linear contact normal/shear damping	${\mathit{\beta }}_{\mathrm{n}}/{\mathit{\beta }}_{\mathrm{s}}$		0.7/0.0
Surface gap	$g_{\mathrm{S}}$	(m)	${\mathrm{2.1}}^{-\mathrm{2}}\times R_{\min }(={\mathrm{5.04}}^{-\mathrm{3}}$)
Gravitational acceleration	$G$	(m s${}^{-\mathrm{2}}$)	9.81

Comparison of cavity generation methods

Several methods have been tested in order to determine the optimal void installation procedure for reasonable simulation time and realistic surface displacement curves. These are instantaneous material removal by particle deletion (M1), incremental material removal by particle deletion (M2), whole cavity particle radii reduction (M3) and incremental particle radii reduction (M4). The radius ($r=\mathrm{5}$ m) and centre depth ($h=\mathrm{35}$ m) of the circular material removal zone was chosen to match the expected sizes of cavities at the area of application. In M1, particles inside the cavity are instantaneously removed, while M2 allows for 15 steps of incremental particle deletion. For the other two methods, parallel bonds are first removed and then we use particle radii reduction in 50 steps to 7.7 % of the original size, each step meaning 5 % reduction. The difference between both is again a complete (M3) vs. incremental (M4) approach. In all methods, the assembly is cycled to SR${}_{\max }=\mathrm{1}e-\mathrm{6}$.

M3 and M4 show similar results for the horizontal displacement $U_x$ but a slightly lower vertical displacement $U_y$ compared with M1 and M2 (Fig. B1). A crucial finding is that M3 and M4 reveal a longer calculation time by 1–2 orders of magnitude. As a result of this test, we consider methods M1 and M2 as generally suitable to simulate a realistic material removal under acceptable calculation time. For the following model verification tests, we rely on method M1 as the simplest option to implement a cavity, and M2 will serve for the final sinkhole models.

Displacement plots for different void installation methods. (a) Horizontal and (b) vertical model setup with size $\mathrm{400}\times \mathrm{400}$ m using a coarse particle distribution (mean radius 0.74 m). Indication of needed simulation time is given in the legend. M1–M4 refer to the installation methods described in the text.

[Figure omitted. See PDF]

Convergence tests on model dimensions and resolution

We performed model resolution tests to determine the optimal size for the mechanical problem of a shallow cavity in a bonded rock assembly. The cavity is installed by instantaneous (quasi-static) particle removal (M1 as shown in Fig. 4a). We varied the width $W$ and height $H$ of the model box from 100 to 800 m while keeping the particle radii constant at 0.74 m for a cavity installation in 35 m depth with a radius of 5 m and track the horizontal and vertical surface displacement.

In Fig. B2, we see the horizontal and vertical displacement curves for all model dimensions. Boundary effects in such a setting close to the free surface make the judgement of the optimal size demanding, but the expected behaviour for the vertical displacement is a subsidence roughly $\mathrm{9}/\mathrm{10}$ and an uplift roughly $\mathrm{1}/\mathrm{10}$ of the total vertical displacement (see model benchmarking in Sect. 2.2).

Convergence test results for model assembly dimensions: cavity of radius of 5 m and depth of 35 m created in each case by method M1. Mean particle radius is 0.74 m. Left and right columns show horizontal and vertical displacement profiles, respectively. Each plot shows results for varying model height ($H$) from 100 to 800 m for a given width ($W$): (a) $W=\mathrm{100}$ m; (b) $W=\mathrm{200}$ m; (c) $W=\mathrm{400}$ m; (d) $W=\mathrm{800}$ m. A convergence is observed for larger model dimensions and a minimum height of 400 m is favoured. Symmetric boundaries ($\mathrm{400}\times \mathrm{400}$ m) give the most stable results.

[Figure omitted. See PDF]

We observe the most stable results for symmetric model dimensions and define the optimal model size to height ($H$) $\times$ width ($W$) $=\mathrm{400}\times \mathrm{400}$ m to account for later possible growth of such a void space. In relation to the cavity size, this means the optimal model is 40 times the cavity diameter. In another expression, the optimal model dimension $/$ cavity depth ratio is 10, a typical value in engineering problems (see Sect. B3). In detail, asymmetric or small model sizes lead to instable results with tails not reaching the expected zero line.

The influence of the particle radii on the displacement curves is shown in Fig. B3 for the above-determined favourable model dimensions. A convergence is observed for particles with mean radius around 0.32 m. Model dimensions of $H\times W=\mathrm{400}\times \mathrm{400}$ m with a mean particle radius of $\stackrel{\mathrm{‾}}{R}=\mathrm{0.32}$ m are thus the optimum parameters to account for converging results, model boundary effects and minimising simulation times.

Results of convergence tests of for influence of particle size: horizontal (a) and vertical displacement (b) profiles for method M1 with in a $\mathrm{400}\times \mathrm{400}$ m sized box. We observe a convergence of the displacement curves for mean radii around 0.32 m, but for decreasing particle sizes, we observe further diminishing of the amplitudes.

[Figure omitted. See PDF]

Detailed of continuum-based solutions for displacements around a gravitationally loaded cavity in 2-D

The first analytical solution used, the Kirsch solution, a classical solution for simple excavation shapes, does not include the free-surface effect, and the mathematical details are depicted, e.g. in Brady and Brown (2006). The radial and tangential displacements at a point $a=a(a,\mathit{\theta }$) at the surface for an average vertical stress $P$, the horizontal stress $K\times P$ and the shear modulus $G$ are $$\begin{array}{cc}{\displaystyle }& {\displaystyle }{u}_r\left(a\mathit{\theta }\right)={\displaystyle \frac{P{r}^{\mathrm{2}}}{\mathrm{4}G\left|a\right|}}\left[\left(\mathrm{1}+K\right)-\left(\mathrm{1}-K\right)\right.\\ & {\displaystyle }& {\displaystyle }\left..\times \left(\mathrm{4}\times \left(\mathrm{1}-\mathit{\upsilon }\right)-{\displaystyle \frac{r^{\mathrm{2}}}{\left|a\right|}}\right)\cos \mathrm{2}\mathrm{\Theta }\right]{\displaystyle }& {\displaystyle }{u}_t\left(a\mathit{\theta }\right)={\displaystyle \frac{P{r}^{\mathrm{2}}}{\mathrm{4}G\left|a\right|}}\left[\left(\mathrm{1}-K\right)\times \left(\mathrm{2}\times \left(\mathrm{1}-\mathrm{2}\times \mathit{\upsilon }\right)\right.\right.\\ & {\displaystyle }& {\displaystyle }\left..\left..+{\displaystyle \frac{r^{\mathrm{2}}}{\left|a\right|}}\right)\sin \mathrm{2}\mathrm{\Theta }.\right]\end{array}$$ With translation into Cartesian coordinates, this yields the surface displacements: $$\begin{array}{}& {\displaystyle }& {\displaystyle }{u}_{xx}\left(y=\mathrm{0}\right)=u_r\cos \mathrm{2}\mathrm{\Theta }-u_t\sin \mathrm{2}\mathit{\theta }& {\displaystyle }& {\displaystyle }{u}_{yy}\left(y=\mathrm{0}\right)=u_r\sin \mathrm{2}\mathit{\theta }+u_t\cos \mathrm{2}\mathrm{\Theta }.\end{array}$$ The second analytical solution used is from Verruijt and Booker (2009) and includes the free-surface effect and is based on the analytical solution of stresses from Mindlin (1940). Verruijt and Booker (2009) added the displacement calculation to the original 2-D Mindlin solution. It is determined via the complex variable method (Muskhelishvili, 1953) and consists of three partial solutions. The second and third partial solutions are relevant for displacement calculation. The second is based on Melan's solution for a concentrated vertical force in a semi-infinite medium and the third involves a balance of the stresses at the cavity boundary. The reader is referred to Verruijt and Booker (2000, 2009) for mathematical details.

The equation for normal displacements as derived by the second solution for an elastic half space ($x$, $y$) under the action of normal line surface load $P$ (Melan's solution) is (Davis and Selvadurai, 1996; Jaeger et al., 2007) $$u_{yy}\left(xy\right)={\displaystyle \frac{\left(\mathrm{1}-\mathit{\upsilon }\right)P}{G\mathit{\pi }}}{\left[\left[\ln (x+a)\right]\right]}_{x=\mathrm{0}}^{x=\mathrm{\infty }},$$ with $G$ the Lamé parameter (shear modulus), $\mathit{\upsilon }$ the Poisson ratio and $a$ the distance to the point of interest.

As known well from linear elastic material theory (Muskhelishvili, 2013; Timoshenko and Goodier, 1973), the integration of the stress formulae is such that a setting of a loaded material (Flamant's problem), which is similar to material removal in the underground, leads to the logarithmic term in the equation above. This leads to infinite vertical displacements along the $x$ surface and a singularity at the centre point ($xy=\mathrm{0}$).

As a workaround for calculation of finite displacements around the cavity, Verruijt and Booker (2009) defined a value $d$ where displacements are set to zero, $u_{yy}(y=d,x=\mathrm{0})=\mathrm{0}$, a so-called fixed point at depth. This constant $d$ can be arbitrarily defined; in engineering, it is usually set to 10 times the depth of the cavity ($d=\mathrm{10}\times h$).

Thus, displacements are considered as not physically realistic in the far field of a load (or cavity), but relative displacement differences are (cf. Davis and Selvadurai, 1996 and Verruijt and Booker, 2009). For the above-stated problem, the relative vertical displacements $\mathrm{\Delta }{u}_{yy}=u_{yy}\left(x_{\mathrm{1}}\right)-u_{yy}(x_{\mathrm{2}}$) between two points $x_{\mathrm{1}}$ and $x_{\mathrm{2}}$ at the surface are (Davis and Selvadurai, 1996) $$\mathrm{\Delta }{u}_{yy}={\displaystyle \frac{\left(\mathrm{1}-\mathit{\upsilon }\right)P}{G\mathit{\pi }}}\ln {\displaystyle \frac{x_{\mathrm{1}}}{x_{\mathrm{2}}}}.$$ Figure B4 highlights the effect of a variation in Young's modulus and the fixed point depth on the fit between modelled vertical displacements and the analytical Mindlin solution described above. A general finding is that $E$ determines the amplitude of the curve and one can gain even better fits of the DEM results when using a higher elasticity module than determined by the simulated rock tests. Furthermore, setting the $d/h$ value to a more realistic value such as 11.43 which corresponds to a cavity central depth of 35 m and a model height of 400 m, shifts the entire vertical displacement curve. The displacement difference is not affected by this integration constant. Hence, when considering the final “best-fit” solution with a low $d/h=\mathrm{4}$ and high elastic modulus $E$ > 10 GPa, the difficulty in cancelling out the integration constant of the analytical displacement solution leaves a still poor fit of the DEM results in the far field but a reasonable fit in the near field of the installed cavity. We use this approach to determine the near field at the surface as approximately $-\mathrm{8}r\le x\le \mathrm{8}r$ with $r$ the radius of the cavity. In our case, this means that the surface near-field limits are $\pm \mathrm{40}$ m from the centre of the depression.

The effect of Young's modulus $E$ and fixed depth point $d/h$ on the vertical displacement (a, c) and displacement difference (b, d) of the Mindlin (MDL) analytical solution as calculated by Verruijt and Booker (2009). Panels (a) and (b) show the scaling effect of the elastic modulus which affects both $U_y$ and $\mathrm{\Delta }{U}_y$. Panels (c) and (d) show the effect of $d/h$ for $E=\mathrm{10}$ GPa which shifts the $U_y$ curve but has no effect on $\mathrm{\Delta }{U}_y$.

[Figure omitted. See PDF]

Details on Mohr–Coulomb and Hoek–Brown rock test analysis

The bulk behaviour of particle assemblies emerges from the interaction of the particle according to the mechanical rules imposed at the contact and bond scale. Therefore, and unlike for continuum-based approaches, the bulk behaviour in DEM models must be calibrated by simulated rock or soil mechanics tests (Potyondy and Cundall, 2004). Here, biaxial compression and tension tests are used to determine the bulk elastic properties of the medium, i.e. the Poisson ratio $\mathit{\nu }$ and Young's elastic modulus $E$. By fitting of the peak stress data upon failure in such tests to, e.g. Mohr–Coulomb or Hoek–Brown failure envelopes, one can also determine bulk strength properties (tensile strength $T$, unconfined compressive strength UCS, coefficient of internal friction $\mathit{\varphi }$).

Differential stress vs. vertical strain for CC and DT tests for a confining pressure of 0.1 MPa: (a) lacustrine mud; (b) alluvium sediments; (c) holocene salt rock. The dashed line indicates elastic limit which was used to determine elastic parameters displayed above the graphs. Red dots mark the peak stresses.

[Figure omitted. See PDF]

A typical stress vs. strain curve contains three parts: (1) a non-linear or linear elastic behaviour, (2) a non-linear yielding behaviour as cracks appear in the material and (3) a non-linear post-peak behaviour after material failure. The peak of the stress–strain curve defines the maximum and minimum principal stresses (${\mathit{\sigma }}_{\mathrm{1}}$, ${\mathit{\sigma }}_{\mathrm{3}}$) at failure. For the compression test, the axial stress is the maximum compressive stress ${\mathit{\sigma }}_{\mathrm{1}}$ (most negative value in the convention used here) and the transversal stress is the minimum compressive stress ${\mathit{\sigma }}_{\mathrm{3}}$ (least negative). For the tension test, it is vice versa: the transversal stress is the maximum tensile stress ${\mathit{\sigma }}_{\mathrm{1}}$ (most positive) and the axial stress is the minimum tensile stress ${\mathit{\sigma }}_{\mathrm{3}}$ (least positive).

The mean peak stresses can be determined for each confining pressure and plotted against each other. In a linear (Mohr–Coulomb) fit of ${\mathit{\sigma }}_{\mathrm{1}}({\mathit{\sigma }}_{\mathrm{3}}$), the UCS is determined by the intercept at ${\mathit{\sigma }}_{\mathrm{3}}=\mathrm{0}$ and the unconfined tensile strength ($T$) by the intercept at ${\mathit{\sigma }}_{\mathrm{1}}=\mathrm{0}.$ The slope $q={\tan }^{\mathrm{2}}\mathit{\beta }$ can be used to fit the Mohr failure envelope as shown in Fig. 7: $${\mathit{\sigma }}_{\mathrm{1}}=C+{\mathit{\sigma }}_{\mathrm{3}}{\tan }^{\mathrm{2}}\mathit{\beta },$$ with $C=\mathrm{UCS}=\mathrm{2}{c}_{\mathrm{0}}\tan \mathit{\beta }$ and $\mathit{\beta }=\mathrm{45}{}^{\circ }+\mathit{\varphi }/\mathrm{2}$. For a Hoek–Brown fit in a ${\mathit{\sigma }}_{\mathrm{1}}({\mathit{\sigma }}_{\mathrm{3}}$) plot, a function of the following form is used: $${\mathit{\sigma }}_{\mathrm{1}}={\mathit{\sigma }}_{\mathrm{3}}+\sqrt{m{\mathit{\sigma }}_{\mathrm{0}}{\mathit{\sigma }}_{\mathrm{3}}+s{{\mathit{\sigma }}_{\mathrm{0}}}^{\mathrm{2}}},$$ with $m$ and $s$ as the empirical rock parameters. For the assumption of intact rock, $s=\mathrm{1}$, ${\mathit{\sigma }}_{\mathrm{0}}=C$, the UCS and $T\sim \frac{C}{m}$. Hence, the fit parameters $m$ and $C$ are used to derive the strength properties of the tested materials. Figure B5 provides exemplary stress vs. strain plots at a confining pressure of $-\mathrm{0.1}$ MPa for all tested materials.

Technical details of implementation of cavity growth and sinkhole collapse in Dead Sea materials in PFC2D v5.0

A PFC- and Python-based code to simulate sinkhole formation

A graphical description of the implemented Python/PFC2D-Fish sinkhole modelling code is depicted in Fig. B6. Here, Fish code parts are marked in yellow and Python code in grey colour. A typical sinkhole simulation follows the following scheme:

• Model dimensions, particle parameters and a function $f\left(i\right)$ for the material removal is defined at the beginning of each set. An unbonded assembly of particles with a fixed porosity of 0.2 is generated at once for the whole assembly at the initial void space growth round ($i=\mathrm{0}$, no material removal yet).

• Similar to the material generation procedure for the model verification material (see Appendix A2), we settle and bond the assembly with a parallel bond model according to the desired material properties. It has to be noted that for low-strength material a bond-reinstallation procedure has been applied, a so-called annealing; i.e. failed bonds can be re-established by contact with other particles of the same material, accounting for, e.g. cohesive mud behaviour. For the other materials, a failed PB is not activated again. We then install the desired tracking functions (measurement circles, marker particles, histories) and group the initial void spaces defined in the model control file.

• This material removal loop acts on each defined cavity growth round $i$. If the area of the particles in the void space zones matches the definition by function $f\left(i\right)$, the loop is broken and important tracked parameters are recorded.

• Step 3 is repeating with increasing material removal round $i$, and after each, the desired tracked results are output via Python code. When a predefined maximum void space growth is reached, the model is finished and a new random assembly starts at step 1.

Graphical description of the PFC2D-based sinkhole modelling code. Yellow colours indicate PFC-Fish language-based code; grey colours are Python control connections. Solid arrows indicate time step cycling. Each model set consists of $n_{\mathrm{rand}}$ random assemblies of particles to account for statistical variation of the DEM sinkhole collapse via arbitrary material removal function in single voids.

[Figure omitted. See PDF]

To avoid another degree of freedom in the calibration of micro- vs. macroproperties, the initial porosity only changes due to the compression by the gravity settling scheme. We have refrained from using either post-settling particle removal to adjust the porosities to specific values or layer-wise gravity deposition with different porosities because of the high amount of calculation time needed.

The Fish material removal core loop (no. 3 in Fig. B6) provides the technical implementation of a quasi-static void space growth. A simple law between the particle area $A_i$ that is supposed to be removed during the void space growth round $i$ and the initial area $A_{\mathrm{0}}$ has been chosen with arbitrarily definable function $f(i$): $$A_i=f\left(i\right)A_{\mathrm{0}}.$$ The void space area is defined by a major and a minor axis. This enables both semi-elliptical, elliptical and circular void space growth. For the results presented in this paper, a slow, constant void space growth $f\left(i\right)={\mathrm{1.0}}^i$ with $A_{\mathrm{0}}=\mathrm{16.3}$ m${}^{\mathrm{2}}$ and a linear eccentricity of $e=\mathrm{2.64}$ was chosen. This avoids the triggering of dynamic effects if too many particles are deleted at once. Other options may include a doubling void space each round ($f\left(i\right)={\mathrm{2.0}}^{i-\mathrm{1}}$ or an exponential increase $f\left(i\right)=e^{(i-\mathrm{1})}$ for $i\ge \mathrm{1}$. For this purpose, a computationally rather cost-intensive static equilibrium procedure is available in PFC2D v5, which sets the bond strengths high before particle deletion, cycles to a stable limit after particle deletion and then resets the bond strengths to the original value.

The pure runtime for a full simulation of an alluvium on mud setup on an Intel Xeon 3.7 GHz processor with 64 GB RAM needs roughly 2 weeks for one particle assembly without tracking geophysical parameters. The tracking would increase the runtime by a factor of $\sim \mathrm{1.5}$. A possible improvement in future will be the introduction of focus regions with an increasing particle radius with distance from the centre of the model.

Details on the implemented parameter tracking

A tracking of pre-, syn- and post-collapse geodetic and geophysical parameters has been implemented in the modelling code (no. 4 in Fig. B6). The technical details are listed as follows.

Porosity, stress and strain rate are recorded using the distribution of so-called measurement circles of area $A^{\mathrm{m}}$ throughout the model domain (Itasca Cooperation Group, 2014; Potyondy and Cundall, 2004).

Porosity is calculated via $n=\frac{V_{\mathrm{void}}}{V_{\mathrm{m}}}=\mathrm{1}-\frac{V_{\mathrm{mat}}}{V_{\mathrm{m}}}$ with $V_{\mathrm{void}}$ as the volume of the void and $V_{\mathrm{mat}}$ the approximated volume of the particles of amount $N_b$ in the measurement circle.

The average stress tensor is calculated in static conditions via $\stackrel{\mathrm{‾}}{\mathit{\sigma }}=-\frac{\mathrm{1}}{V_{\mathrm{m}}}\sum _{N_c}{F}^{\left(\mathrm{c}\right)}\otimes L^{\left(\mathrm{c}\right)}$ where $\otimes$ is the dyadic product of two tensors, $N_{\mathrm{c}}$ is the number of contacts, $F^{\left(\mathrm{c}\right)}$ is the contact force vector and $L^{\left(\mathrm{c}\right)}$ the branch vector that joins the centroids of two entities. From this, the maximum compression principal stress ${\mathit{\sigma }}_{\mathrm{1}}$, the minimum compression principal stress ${\mathit{\sigma }}_{\mathrm{3}}$ and the maximum shear stress $T_{\max }=\frac{({\mathit{\sigma }}_{\mathrm{3}}-{\mathit{\sigma }}_{\mathrm{1}})}{\mathrm{2}}$ are calculated, which is always positive in the convention used here, where compression is negative.

The strain rate tensor $\dot{e}$ (velocity gradient tensor) is calculated via a least-squares best-fit approach of the predicted vs. the measured relative velocities ${\stackrel{\mathrm{̃}}{V}}_i^{\left(p\right)}=V_i^{\left(p\right)}-{\stackrel{\mathrm{‾}}{V}}_i$ between each two entities $p$ during a time step $i$, with $V_i^{\left(p\right)}$ as the translational velocity and ${\stackrel{\mathrm{‾}}{V}}_i=\frac{{\sum }_{N_p}{V}_i^{\left(p\right)}}{N_p}$ as the mean velocity in the circle area.

The strain tensor in the measurement region is then calculated by multiplying strain rate components with the simulation time step and summing over the desired period.

Alternatively, strain is calculated via simulated extensometers. For these, pairs of particles which lie either horizontally or vertically next to each other are defined. By registering the displacement of each particle, a pairwise calculation of the horizontal and vertical strains is achieved at low computational cost in comparison to the measurement circle distribution (Itasca Cooperation Group, 2014).

A full set of metadata is available upon request. For photogrammetric surveys, raw images, DSMs and orthophotos are available upon consultation with the authors. For DEM models, data and results are available upon request.

The supplement related to this article is available online at: https://doi.org/10.5194/se-9-1341-2018-supplement.

DAH led the production of figures, data analysis and numerical modelling. DAH and EPH led the writing of the manuscript. DAH, EPH, AT, MPJS and SE developed the calibration of the materials and discussed technical issues. All authors reviewed and commented on the manuscript, and they contributed to discussions of the data and interpretation.

The authors declare that they have no conflict of interest.

This article is part of the special issue “Environmental changes and hazards in the Dead Sea region (NHESS/ACP/HESS/SE inter-journal SI)”. It is not associated with a conference.

Acknowledgements

We would like to thank the two reviewers and the editor for the fruitful discussion. We very gratefully acknowledge fieldwork support from Ali Sawarieh, Hussam Alrshdan and their colleagues at the Ministry of Energy and Mineral Resources of the Hashemite Kingdom of Jordan. Acknowledgements go to Emad Talafheh and Zayad Mansour from the Arab Potash Company, as well as to Michael Ezersky of the Geophysical Institute of Israel, for their support regarding geotechnical parameters. Also, we thank Arnold Verruijt for the comments regarding the analytical solution, as well as Daniel Woodell and Mehdi Nikhoo for the technical advice of the DEM and BEM modelling, respectively, and Marc Elmouttie for facilitating the research stay at CSIRO (Brisbane). Special thanks go to Damien Closson and Thomas R. Walter for continuous support and material provision. Particular thanks go to Itasca Consulting Group for providing the license of PFC2D v5.0 in the framework of the Itasca Education Partnership program and the German Academic Exchange Service (DAAD) for a short-term doctorate research grant. Last, but not least, thanks go to the projects DESERVE and SIMULTAN and the involved colleagues for their kind support and funding opportunities. The article processing charges for this open-access publication were covered by a Research Centre of the Helmholtz Association. Edited by: Amotz Agnon Reviewed by: Renaud Toussaint and one anonymous referee