Dynamics of hydrological and geomorphological processes in evaporite karst at the eastern Dead Sea – a multidisciplinary study

Karst groundwater systems are characterized by the presence of multiple porosity types. Of these, subsurface conduits that facilitate concentrated, heterogeneous flow are challenging to resolve geologically and geophysically. This is especially the case in evaporite karst systems, such as those present on the shores of the Dead Sea, where rapid geomorphological changes are linked to a fall in base level by over 35 m since 1967. Here we combine field observations, remote-sensing analysis, and multiple geophysical surveying methods (shear wave reflection seismics, electrical resistivity tomography, ERT, self-potential, SP, and ground-penetrating radar, GPR) to investigate the nature of subsurface groundwater flow and its interaction with hypersaline Dead Sea water on the rapidly retreating eastern shoreline, near Ghor Al-Haditha in Jordan. Remote-sensing data highlight links between the evolution of surface stream channels fed by groundwater springs and the development of surface subsidence patterns over a 25-year period. ERT and SP data from the head of one groundwater-fed channel adjacent to the former lakeshore show anomalies that point to concentrated, multidirectional water flow in conduits located in the shallow subsurface (< 25 m depth). ERT surveys further inland show anomalies that are coincident with the axis of a major depression and that we interpret as representing subsurface water flow. Low-frequency GPR surveys reveal the limit between unsaturated and saturated zones (< 30 m depth) surrounding the main depression area. Shear wave seismic reflection data nearly 1 km further inland reveal buried paleochannels within alluvial fan deposits, which we interpret as pathways for groundwater flow from the main wadi in the area towards the springs feeding the surface streams. Finally, simulations of density-driven flow of hypersaline and undersaturated groundwaters in response to base-level fall perform realistically if they include the generation of karst conduits near the shoreline. The combined approaches lead to a refined conceptual model of the hydrological and geomorphological processes developed at this part of the Dead Sea, whereby matrix flow through the superficial aquifer inland transitions to conduit flow nearer the shore where evaporite deposits are encountered. These conduits play a key role in the development of springs, stream channels and subsidence across the study area.

1 Introduction

Karst landscapes result from the dissolution by water of rocks or semi-consolidated sediments on and below the Earth's surface (Ford and Williams, 2007). Typically, weakly acidic meteoric water percolates into the subsurface, where it dissolves soluble material and enhances its porosity beyond initial matrix porosity and/or fracture porosity (Price, 2013; Hartmann et al., 2014; Goldscheider, 2015). Consequently, a characteristically well-evolved network of subsurface voids and conduits in karst groundwater systems enhances the hydraulic conductivity of the host materials to extremes rarely seen in non-karstified aquifers (Kaufmann and Braun, 2000; Hartmann et al., 2014).

The scale and geometry of karst conduit networks are variable, and they depend upon factors including the lithological characteristics of the host rock, the regional geological and climatic setting and the nature of hydrological recharge (Ford and Williams, 2007; Gutiérrez et al., 2014; Parise et al., 2018). In limestone areas of high annual precipitation, for example, surface dissolution leads to the formation of enclosed depressions (“solution dolines”), which channel surface water into the underground (Waltham et al., 2005; Sauro, 2012). Fractures formed along faults, joints and bedding planes within the host rock can be enlarged by dissolution to create branching networks of caves which are large enough to be entered by speleologists (Palmer, 1991, 2007, 2012; Ford and Williams, 2007). As a result of their global prevalence and ready access, karst systems formed in fractured limestone constitute the bulk of our present understanding of karst drainage systems.

In contrast, karst systems formed by dissolution of evaporites are less common and less studied. Evaporite deposits are commonly poorly bedded and may form constituent parts of marine or lacustrine sedimentary deposits that are semi-consolidated to fully consolidated (Warren, 2006). Due to the extreme solubility of evaporite minerals such as halite, they are only able to form in arid environments (Frumkin, 2013), such that surface recharge is limited as compared to typical limestone karst. Moreover, karst systems in young evaporites are characterized by a more dynamic, changing flow system due to the instability of the flow tubes (Ford and Williams, 2007). Baseflow in these conduits may be of very low discharge, and surface flow may only occur rarely in ephemeral wadis (dry river valleys), where water often flows beneath the surface (Goode et al., 2013; Price, 2013; Salameh et al., 2018).

Figure 1 Overview of the Ghor Al-Haditha field site. (a) Topographic map showing the location of the study area and field site on the eastern Dead Sea shore, Jordan. Elevation data are derived from the Shuttle Radar Topography Mission (Farr et al., 2007). The inset shows the catchments of the three major wadis flowing to the Dead Sea shoreline in the study area: Wadi Mutayl, Wadi Ibn Hamad, and Wadi Al-Mazra'a. (b) Pleiades 2018 satellite image of the study area showing the studied groundwater-fed stream channels (blue) mapped from satellite imagery and close-range photogrammetric surveys. Label CM- refers to meandering-type channels, and CS- refers to straight or braided channels. Geophysical survey lines or areas, as analysed in this work, are marked. Seismic lines SL1 and SL2 refer to S-wave reflection profiles 1 and 2 acquired in 2014 (Polom et al., 2018). The two boreholes of El-Isa et al. (1995) are labelled “BH1” and “BH2”. The area affected by sinkhole formation from 1985 to present is represented by the dashed outline with shaded infill. The numbered blue and white dots represent the approximate locations of groundwater springs sampled during a previous field campaign in the study area in 1999, from Sawarieh et al. (2000), and discussed further in Sect. 2.

One of the most rapidly developing evaporite karst systems on the globe is found on the shores of the Dead Sea (Fig. 1a). This hypersaline terminal lake, fed primarily by the Jordan River, lies within the ∼ 150 km long and ∼ 8–15 km wide Dead Sea basin (Garfunkel and Ben-Avraham, 1996). This basin has subsided rapidly from the late Pliocene to present (Ten Brink and Flores, 2012) due to motion along the left-lateral Dead Sea Transform fault (DSTF) system. During this time, several paleolakes of varying size and duration (Bartov, 2002; Torfstein et al., 2009) have existed within the basin. Since the 1960s, the modern Dead Sea level has declined from −395 to −434 $\mathrm{m}\mathrm{m}.\mathrm{s}.\mathrm{l}.$ (1967–2020; ISRAMAR, 2020), i.e. by 40 m as of 2020. The lake level fell by 0.5 m yr⁻¹ in the 1970s and by 1.1 m yr⁻¹ in the last decade. This has led to a dynamic reaction of the hydrogeological system and landscape around the lake shore (Kiro et al., 2008). New springs, stream channels, depressions, landslides and sinkholes have formed since the 1980s. These have developed both within alluvial fan sediments deposited by flash floods in wadis terminating at the lake and within mud and evaporite deposits of the former lake bed that have been revealed by the retreat of the shoreline (see e.g. Abelson et al., 2006, 2017; Closson et al., 2007; Parise et al., 2015; Yechieli et al., 2016; Al-Halbouni et al., 2017; Watson et al., 2019, and references therein). The hazards posed by erosion and subsidence result in severe consequences for tourism, agriculture and infrastructure at the Dead Sea (Arkin and Gilat, 2000; Closson and Abou Karaki, 2009; Fiaschi et al., 2018; Abou Karaki et al., 2019).

There are several hypotheses regarding the nature of the subsurface hydrogeology at the Dead Sea and how this relates to the development of rapid erosion and subsidence around the shoreline. There is a consensus that the falling Dead Sea level increases the hydrological gradient adjacent to the lake. This promotes deeper incision of existing stream channels around the lake shore, adjustments in the planform morphology of existing channels (e.g. increased sinuosity), and incision of new channels on the exposed lake bed (Ben Moshe et al., 2008; Bowman et al., 2007, 2010; Dente et al., 2017, 2019; Vachtman and Laronne, 2013). Several studies contend also that the base-level fall results in lateral migration of the freshwater–saltwater interfaces in the subsurface, which enables relatively undersaturated groundwater to invade evaporite deposits previously in equilibrium with Dead Sea brine, thus driving karstification of those deposits and subsidence of the ground surface (Salameh and El-Naser, 2000; Yechieli, 2000; Yechieli et al., 2006). Other studies have argued that karstification is driven primarily (or additionally) by preferential upward flow of fresher groundwater into the evaporite-rich deposits via hydraulically conductive regional tectonic faults (Abelson et al., 2003; Closson et al., 2005; Shalev et al., 2006; Closson and Abou Karaki, 2009; Charrach, 2018). Two main approaches have been used to simulate mathematically such hydrogeological scenarios in or near the DS rift valley: (1) a “sharp-interface approximation” to the transition between hypersaline and fresh groundwaters (Salameh and El-Naser, 2000; Yechieli, 2000; Kafri et al., 2007) and (2) density-driven flow simulations that account for mixing of waters of contrasting salinities and densities (e.g. Shalev et al., 2006; Strey, 2014). Both approaches yield some success within the existing hydrogeological constraints.

There are also varying views of the form and make-up of the subsurface evaporite deposits and varying emphasis on the mechanisms of subsurface erosion that lead to surface subsidence. Many studies invoke chemical erosion (dissolution) of a massive salt layer of up to 15 m thickness and lying at depths of 20–50 m. This layer is proposed to be dissolved by diffuse or fracture-controlled groundwater flow, which results in large intra-salt cavities that then collapse to produce sinkholes (Shalev et al., 2006; Ezersky and Frumkin, 2013; Yechieli et al., 2016; Abelson et al., 2017). Other studies emphasize physical erosion (“piping” or “subrosion”) of weakly consolidated alluvial sand and gravel deposits due to focussed, turbulent groundwater flow in such materials at depths of 5–20 m (Taqieddin et al., 2000; Sawarieh et al., 2000; Arkin and Gilat, 2000). Such studies further propose a combined physical and chemical erosion of thinly interbedded lacustrine salt, marl and clay deposits at the depths of 1–40 m, whereby salt dissolution and clay weakening by relatively fresh groundwater flow generate subsurface cavities and conduits (Arkin and Gilat, 2000; Al-Halbouni et al., 2017; Polom et al., 2018). It is emphasized here that such varying views are not mutually exclusive (cf. Watson et al., 2019). Indeed, some recent studies have emphasized linkages between the surface water flow, groundwater flow in a shallow conduit system and dissolution of deeper-seated salt – particularly during flash-flood events – as an important factor in sinkhole formation (Avni et al., 2016; Arav et al., 2020).

In consensus, however, all mentioned processes lead to mechanical weakening of the underground materials and to subsidence, which is expressed as sinkholes of 1–70 m width and up to 15 m depth and as uvala-like depressions of several hundred metres in width and up to 12 m in depth (Watson et al. 2019). The style of subsidence depends strongly on the mechanical properties of the materials involved (Al-Halbouni et al., 2018, 2019). Both cover-collapse and cover-sag sinkholes occur, according to the classification of Gutiérrez et al. (2014) and Parise (2019).

The use of indirect (non-invasive) methods to investigate subsurface hydrogeology in karst environments presents unique methodological challenges due to the high degree of heterogeneity and discontinuity in the subsurface (Palmer, 2007; Goldscheider and Drew, 2007; Goldscheider, 2015). An overview of indirect methods applied to the Dead Sea sinkhole phenomenon is given by Ezersky et al. (2017). Despite these challenges, near-surface geophysical methods have shown potential for cost-effective investigation of subsurface erosion in evaporite karst environments (see e.g. Wust-Bloch and Joswig, 2006; Jardani et al., 2007; Gutiérrez et al., 2011; Krawczyk et al., 2012; Malehmir et al., 2016; Giampaolo et al., 2016; Fabregat et al., 2017; Abelson et al., 2018; Muzirafuti et al., 2020) and determination of the actual configuration of the fresh–saline water interface (e.g. Kruse et al., 2000; Kafri and Goldman, 2005; Kafri et al., 2008; Meqbel et al., 2013). This is especially the case when geophysical surveys are supported by borehole investigations and other lines of evidence to aid interpretations. Geophysical investigations combined with borehole and hydrogeochemical data on the western shore of the Dead Sea, for instance, provide strong evidence of dissolution of a massive salt layer in the shallow subsurface (as described above) with cavity formation as imaged by high geoelectric resistivities and for linkage of these processes to sinkhole development in many areas there (e.g. Yechieli et al., 2006; Ezersky, 2008; Ezersky et al., 2011). Such evidence is weaker or more ambivalent on the eastern shore, however (Al-Zoubi et al., 2007; Frumkin et al., 2011; Polom et al., 2018).

In this paper, we investigate the relationship between geomorphological change and subsurface groundwater flow on the Dead Sea's eastern shoreline, at a key part of the Ghor Al-Haditha sinkhole site in Jordan. We combine surface manifestations, remote-sensing analysis, and multiple geophysical surveying methods of horizontally polarized shear wave (S_H) reflection seismics, electrical resistivity tomography (ERT), self-potential (SP) and ground-penetrating radar (GPR). These are complemented by numerical simulations of the density-driven groundwater flow.

First, we summarize the regional and local geology to aid interpretation of geophysical data and to assign realistic modelling conditions. We then use remote-sensing and surface manifestations to describe the spatio-temporal evolution of a set of stream channels fed by groundwater to show a major shift in the local hydrology associated with numerous sinkhole collapses, particularly at the head of a late-developed major channel formed in the exposed lacustrine deposits at the Dead Sea's eastern shoreline. We then present new ERT and SP datasets gathered at the head of this major channel to image the subsurface flow paths, as well as new ERT data that complement previous shear wave seismic survey data (Polom et al., 2018). We also present new shear wave seismic data further inland, near to the outlet of the main wadi terminating in the area. We provide constraint on the groundwater table by GPR surveys. On the basis of all the above, we hypothesize that flow of subsurface groundwater undergoes a transition from a weakly focused matrix flow near the main wadis to strongly focussed conduit flow near to the shore before feeding the major channel. To test this, we present 2D density-driven groundwater flow modelling that suggests that such conduits are a necessary component of the groundwater system to generate realistic groundwater head and electrical conductivity (EC) conditions in the case of a regressing Dead Sea level. This multidisciplinary analysis leads us to present a refined picture of the groundwater system at Ghor Al-Haditha with improved insight into the links between surface subsidence and erosion and groundwater flow in the case of the regressing Dead Sea.

Figure 2 Field photos of the Ghor Al-Haditha study area. (a) View from the Moab mountains near Wadi Abu-Darat onto the former Dead Sea lake bed. The dashed white line marks the boundary between alluvium and lacustrine deposits, and it is approximately equivalent to the shoreline position in 1960. (b) Vice versa view from the DS lake bed towards the Moab mountains, with a section of a meandering stream channel and associated deformed ground. (c) Area around canyon CM5 with the destroyed Numeira mud factory. (d) Zoom-in of lacustrine deposits. (e) Zoom-in of alluvial fan deposits.

2 Site of investigation and previous geophysical and hydrogeological studies there

The Ghor Al-Haditha study area is situated on the eastern shore of the Dead Sea (Fig. 1). The climate is semi-arid to arid, with average annual precipitation of 70–100 mm (El-Isa et al., 1995; Salameh et al., 2018). Three major wadi systems, Wadi Ibn Hamad, Wadi Mutayl and Wadi Al-Mazra'a (or Wadi Al-Karak) terminating in the area, have deposited sequences of semi-consolidated to unconsolidated alluvial fan deposits at the coastline (Figs. 1b and 2a). The drainage basins of the wadis are delineated in the inset of Fig. 1. The catchment areas of the wadis are 30 km² for Wadi Mutayl, 124 km² for Wadi Ibn Hamad, and 188 km² for Wadi Al-Mazra'a. West and north of the alluvial plain, exposure of the former lake bed by the Dead Sea's recession has formed a “mudflat” or “salt flat”.

Figure 3 Overview of the hydrogeology of the Ghor Al-Haditha study area. (a) Simplified geological map of the study area, partly based on 1:50 000 scale mapping of the Jordanian Ministry of Energy and Mineral Resources (Khalil, 1992), mapping and reports (IOH and APC, 1995) and partly our own work. The stratigraphy generally dips acutely to the south-east while striking to the north-east. Also shown is the right-lateral oblique Siwaqa fault, the inferred position of the Eastern Boundary Fault (downthrowing to the west), and the axis of the Haditha syncline. The present locations of the groundwater spring rising in the Wadi Sir limestones of the Ajlun group, “Ain Maghara” and the spring rising at the contact between alluvium and mudflat deposits at the head of channel CM5 are also shown along with the location of the large uvala-like depression described in Al-Halbouni et al. (2017) and Watson et al. (2019). A 2D cross section of geology along the line A–A` is shown in (b). The red and black star indicates the position of the Numeira mud factory, which has been defunct since 2009 after being destroyed by sinkhole formation at the site. (b) Cross section of vertical exaggeration x5 depicting the subsurface hydrogeology as reported for the early months of 1995 (El-Isa et al., 1995) based on borehole information, refraction seismic surveys and VES surveys. The projected depth to the Ram-Kurnub sandstones is based on borehole data (IOH and APC, 1995). (c) Inset depicting the hypothesized conceptual model of seaward migration of sinkhole populations at the Dead Sea as the fresh–saline interface migrates with the base-level drop (Salameh and El-Naser, 2000; Yechieli, 2000; Abelson et al., 2003, 2017; Watson et al., 2019). “GW” here refers to the inferred level of the groundwater in the subsurface; squiggly arrows show the movement of water in the subsurface.

The lacustrine deposits comprise alternating, light–dark, laminated to thin layers of carbonates (aragonite, calcite), quartz and clay minerals (kaolinite), with centimetre-scale, idiomorphic halite crystals. These are interbedded with localized thin to thick beds of halite and/or gypsum (Fig. 2b and d), the proportion and thicknesses of which near or at the surface increase lakeward and northward (Al-Halbouni et al., 2017; Watson et al., 2019). The alluvium consists generally of poorly sorted, semi-consolidated to unconsolidated sands and gravels (Fig. 2c and e), which are interbedded with minor silts and clays (El-Isa et al., 1995; Taqieddin et al., 2000; Sawarieh et al., 2000; Polom et al., 2018). At the former stable shoreline, alluvium and lacustrine deposits are in direct lateral contact. Comparable deposits have been widely found along the western shore of the Dead Sea and in numerous boreholes (Yechieli et al., 1993, 2002). Data for the Ghor Al-Haditha study area are considerably less comprehensive, consisting of two boreholes, 45 and 51 m deep, drilled on the alluvial plain in January–February of 1995 (El-Isa et al., 1995; BH1 and BH2, Figs. 1b and 3). Borehole logs show alluvium down to the bottom metre, where a layer of silt and greenish clay was detected. The exact ages of the alluvium and lacustrine deposits are unclear; the exposed lacustrine deposits likely correspond to the regional Ze'elim formation of Holocene age (Yechieli et al., 1993), while much of the alluvium may belong to the Lisan formation of Pliocene–Pleistocene age (Khalil, 1992).

There are three principal aquifer units in the area (Fig. 3) whose characteristics are summed up in the following based on works of Akawwi et al. (2009), Goode et al. (2013), IOH and APC (1995), and Khalil (1992). The first is a lower sandstone aquifer comprising the Ram group and Kurnub formation of Cambrian to early Cretaceous ages, respectively, which is not present at the surface in the study area but outcrops along the escarpment of the DSTF just to the north (cf. Fig. 1c, Watson et al., 2019). The second is an upper carbonate aquifer spanning the Ajlun and Belqa groups of late Cretaceous to early Tertiary age. The third is a superficial aquifer in the Lisan and/or Ze'elim formations. Recharge to the Ajlun-Belqa and superficial aquifers is primarily derived from precipitation in the highlands to the east, where the average annual precipitation is ∼ 350 mm. As on the western side of the Dead Sea (Yechieli et al., 2006), it is estimated that some recharge occurs also via leakage from the lower Ram-Kurnub aquifer (cf. p. 5, Sect. 2.3, IOH and APC, 1995). From these regional aquifers, groundwater flows toward the Dead Sea. A significant spring within the Ajlun group limestones, “Ain Maghara”, can be found just to the east of the study area (Fig. 3). The discharge of this spring is around 0.28 m³ s⁻¹ and is thought to be principally derived (∼ 80 %) from baseflow within the Wadi Ibn Hamad (cf. p. 18, Sect. 5, IOH and APC, 1995), which in turn originates from outcrops of the Ram-Kurnub aquifer unit along the wadi bed upstream of Ain Maghara. The remaining (∼ 20 %) flow is thought to derive from upward leakage from the Ram-Kurnub aquifer. Several springs are found in the transition between the alluvium and lacustrine deposits; these feed surface streams that drain into the Dead Sea via numerous channels (Fig. 1b and 2b). Other surface stream channels in the lacustrine deposits lack these groundwater-fed springs and are instead fed by surface water from wadis during flash-flood events. A borehole drilled in the wadi bed just upstream of Ain Maghara confirms subsurface water flow within the wadi bed with a similar chemistry to that of the spring, though the depth and nature of this flow are not recorded. Head elevations in the vicinity of Ain Maghara were reported to be −300 to −350 $\mathrm{m}\mathrm{m}.\mathrm{s}.\mathrm{l}.$ in 1994 (cf. p. 5, Sect. 2.2, IOH and APC, 1995). The depth to the water table in the superficial aquifer at the time of drilling in 1994 was 20.5 m in BH1; no depth was apparently recorded for BH2 (El-Isa et al., 1995). A groundwater well just over 5 km south of the study area, at the Al-Mazra'a pumping centre, indicates that the groundwater head level in the superficial aquifer declined at an average rate of 0.75 m yr⁻¹ from 1960 to 2010, whilst over the same time period the measured electrical conductivity had increased by 59 µS yr⁻¹ (Goode et al., 2013).

The first published report of geophysical investigations in the study area is El-Isa et al. (1995), which includes results of classic refraction seismic and vertical electrical sounding (VES) surveys conducted in 1994. Several anomalies present in the refraction seismic data were inferred to represent buried alluvial channel deposits at 5–10 m depth, just south-west of BH2 (in the centre of their profile 4 and at the southern end of their profile 5; cf. Figs. 4.9 and 4.11 in El-Isa et al., 1995). Ten VES traverses were performed in the study area using an Atlas-Copco SAS-300 system. The combined inversion of survey data across all traverses, calibrated with respect to the hydraulic head measured in the boreholes, enabled El-Isa et al. (1995) to create a layered 2D interpretation of groundwater conditions, consisting of an uppermost “dry” area (resistivities of 120–660 Ωm), a fresh groundwater layer (resistivities of 10–80 Ωm), a brackish zone where groundwater and Dead Sea water mixing occurs (resistivities of 2–8 Ωm) and then saline Dead Sea water (resistivities of < 1 Ωm) and is sketched in Fig. 3b. The water levels in the boreholes were combined with spring elevation data to determine a hydraulic gradient of roughly 30 m km⁻¹ from east to west, following the surface topography.

Table 1 Water geochemistry of samples taken from groundwater springs during the 1999 field campaign of Sawarieh et al. (2000). Spring ID numbering is that used in the 1999 field campaign report; spring locations are shown approximately in Fig. 1b.

A second report, by Sawarieh et al. (2000), includes additional refraction seismic surveys, as well as water chemistry, temperature, pH and EC of springs, wells and sinkholes undertaken in February 1999 (see Fig. 1b for the locations of the springs sampled). The seismic results generally resolved an up to three-layer velocity model, thought to represent differing levels of compaction of alluvium. The hydrology of the study area has changed since the survey: several springs which they sampled in the area of Wadi Mutayl have now dried up, and the sampled locations do not match those of the groundwater springs feeding channels CM1–3 and CM6 (Fig. 1b). Groundwater samples were analysed for anion and cation concentrations and total dissolved solids (TDS) in the laboratory at the Ministry of Energy and Mineral Resources of Jordan (MEMR). The results are summarized in Table 1. In general, the measured TDS increased from east to west: groundwater sampled at Ain Maghara could be defined as between fresh and brackish (EC: ∼ 0.7–0.9 mS cm⁻¹; TDS: ∼ 500–700 ppm), whereas water rising from the springs closer to the shore would be categorized as brackish (EC: 1.5–50 mS cm⁻¹; TDS: 500–35 000 ppm), aside from their sample 1, which had abnormally high EC and TDS, being more of a brine. The $\mathrm{Na}/\mathrm{Cl}$ molar ratios of water samples from these springs were relatively low. The EC readings from the field broadly correlate with the proportions of TDS and Na⁺ and Cl⁻ ions measured in the water samples.

From 2009 to 2010, scientists of MEMR performed geophysical surveys using ERT, GPR, and transient electromagnetic techniques (TEM, Alrshdan, 2012). GPR data acquired with a SIR-10B system near the mud factory were interpreted as very shallow (< 7 m) sinkhole and subsurface water conduits; the same relates to surveys further inland. The used 100 MHz antenna and wet and saline ground conditions strongly restricted the penetration depth of the GPR signal (Alrshdan, 2012). Similar GPR investigation of the shallow subsurface and sinkhole formation was reported by Batayneh et al. (2002). ERT surveys were performed traversing the alluvium–mudflat boundary close to CM1 and the former mud factory site and further inland in farmers' fields parallel to the road along which ERT1 was performed in this study. The survey lines near the mud factory imaged some extremely high-resistivity anomalies (> 10 000 Ωm) thought to represent subsurface cavities. Three ERT profiles were acquired with the IRIS Syscal system by Al-Zoubi et al. (2007), in the sinkhole area NE of Wadi Ibn Hamad. The higher resistivities up to 600 Ωm were interpreted as fractured alluvial material, while low resistivities of less than 20 Ωm were recorded at the bottom of all the profiles. Further ERT, GPR, TEM, magnetic and seismic measurements including multichannel surface wave analysis were performed on the agricultural fields (Camerlynck et al., 2012; Bodet et al., 2010), the seismics part of which has been reported by Ezersky et al. (2013a, b). The two ERT profiles acquired further north in the study area also imaged high-resistivity anomalies (Frumkin et al., 2011), which they interpreted as representing subsurface cavities. All ERT and GPR surveys on the alluvial fan deposits only penetrated to a maximum depth of 25 m, which was insufficient to image the water table in this area.

The most recent geophysical survey of Polom et al. (2018) comprises a summary of shear wave reflection seismic data collected in 2013/14. These data were interpreted as incompatible with a proposed 2–10 m massive salt layer at around 35–40 m depth (Ezersky et al., 2013a) but did reveal zones of high S-wave scattering in the subsurface near the main depression as well as potential buried and refilled channel systems within the alluvial fan architecture and a deeper-lying silt and clay layer at around 60–80 m depth.

3 Material and methods

Our multidisciplinary investigation of surface and subsurface hydrology at Ghor Al-Haditha in Jordan consists of (1) remote sensing based on satellite imagery and aerial photogrammetric image analysis, (2) geophysical methods such as ERT, SP, GPR and S_H-wave reflection seismic surveys, and (3) 2D hydrogeological modelling of density-driven groundwater flow. Detailed theoretical concepts and additional results of the methods are illustrated in Appendix A.

3.1 Remote-sensing and field surveys

The remote-sensing dataset comprises high-resolution optical satellite imagery and aerial photographs that span the years 1967 to 2018. The dataset is identical to that presented in Table 1 of Watson et al. (2019), save for an additional Pleiades optical satellite image acquired on 23 April 2018. Temporally, the satellite and aerial imagery dataset is decadal in resolution from 1967 and annual in resolution from 2004, with spatial resolutions from 2 to 0.3 m pixel⁻¹. The remote-sensing data are complemented by very high-resolution (0.1 m pixel⁻¹) orthophotos and digital surface models derived by photogrammetric processing of optical images obtained through balloon- and drone-based close-range aerial surveys of the study area in 2014, 2015 and 2016 (cf. Table 1, Watson et al., 2019).

The satellite images were pre-processed (orthorectification, pansharpening and georeferencing and co-registration) by using the PCI Geomatica Orthoengine software package, as described in Watson et al. (2019), except for the 2016, 2017 and 2018 Pleiades images which were orthorectified and georeferenced by Airbus. For details of the method for generation of the orthophotos and DSMs from photogrammetric survey data, see Al-Halbouni et al. (2017). After pre-processing, all remote-sensing and photogrammetric datasets were integrated within the Q-GIS Geographical Information Systems (GIS) software package for geographical analysis, which included the manual digitization of fluvial and karst features in order to reconstruct their development through both space and time.

Maps derived from remote-sensing and photogrammetric data were extensively ground-truthed over the course of four field campaigns in October 2014, October 2015, December 2016 and October 2018. The drone and balloon surveys were performed in conjunction with the first three field campaigns. In 2015, we also made a preliminary survey of water bodies in the study area. Temperature and electric conductivity were measured in situ at numerous springs and ponds (within sinkholes) with a Hach HQ40D Portable Multi-meter. These measurements were repeated for selected springs and ponds during the 2018 field campaign.

3.2 Geophysical surveys

3.2.1 Self-potential

The self-potential (SP) method is a passive geophysical method that measures the electric potential between two non-polarizing electrodes connected on the ground surface of the Earth. More information on the physical background of the SP method is given in Appendix A.

Figure 4 Ground-based images of the ERT, SP and seismics equipment. (a) Seismics line 1 with landstreamer geophone system and ELVIS S_H shear wave generator system. (b) ERT line 3 at the head of canyon CM5 with electrode chain ActEle spread. (c) ERT Lippmann 4point light and DGPS Trimble ProXRT-2 equipment. (d) Close-up of stainless-steel electrode, connecting cable and switch box of the multi-electrode ERT system. (e) SP Ag–AgCl electrode in a cotton bag filled with bentonite and watered to improve electric ground coupling. (f) GPR Malå Ramac device with a 50 MHz transmitter and receiver antennas.

Our SP survey was conducted in 2018 (for the location, see Fig. 1b). We used the fixed base configuration (Fig. A1a), with a high impedance (> 50 MΩ) Voltcraft voltmeter and several non-polarizable Ag/AgCl electrodes for the survey (Fig. 4e). The semi-permeable bottom filters of the electrodes are covered by a thick, wet bentonite mass inserted into a closed cotton bag. In this configuration, the electric potential is measured between a fixed base reference electrode and a moving electrode. Adapted to our survey area, we measured parallel profiles to form an array configuration with an inter-electrode spacing of 10 m. The reference electrode, with the convention that the negative pole should be N or E of the positive pole, is buried and watered by using freshwater to improve electrical ground coupling. The moving electrodes are also watered, and after the signal stabilized, we took at least three measurements for each point. More details on the SP array can be found in the results (Sect. 4.2.1).

3.2.2 Electric resistivity tomography

Electric resistivity tomography (ERT) is an active geophysical method within the family of geoelectric methods. Geoelectrics were developed initially as vertical electric sounding (1D) to determine groundwater in the subsurface (Schlumberger, 1920) by injection and recording of an electric field using two pairs of electrodes. More details on the physical background of ERT are given in Appendix A.

Our ERT survey was performed in October 2018 by using a Lippmann 4point light 10 W direct current instrument (Grinat et al., 2010) and stainless-steel electrodes connected to the Lippmann multichannel cable system (Fig. 4b–d). Data were collected by using a multichannel control system implemented in the Geotest v. 2.46 software (Rauen, 2016) with automatic combination of all possible injection (AB) and measuring (MN) pairs of electrodes using a measurement interval frequency of 8 Hz. Different survey configurations for certain field conditions and penetration depths were tested on site prior to the survey. Multichannel geoelectric control systems commonly enable a combination of different general electrode configurations, i.e. Wenner, Schlumberger and several dipole–dipole configurations. The multichannel control system used here was a combination of Wenner and Schlumberger configurations (Fig. A1b) which incorporates the advantages of both general configurations regarding advantages of lateral mapping and deep sounding. In profiles along the roads (ERT1 and ERT2) we used the so-called roll-along technique to provide a combination of long profile distances (0–600 m) with short inter-electrode distances (1–5 m). Repeated measurements have been performed at certain electrode positions in case of non-reliable high-resistance (> 1 KΩ) results and with the support of freshwater to improve electrical ground coupling. Remaining anomalous data points (showing e.g. physically implausible negative resistivities or non-reliable local high resistivities) have been carefully removed from further analysis.

3.2.3 Shear wave reflection seismics

The shear wave reflection seismics method uses the principles of reflection seismology (Sheriff and Geldart, 1995), which is widely applied in hydrocarbon exploration and explained more in Appendix A.

The shear wave seismic surveys reported on here used horizontally polarized S waves (S_H waves) and were undertaken in October 2014 and December 2016. Details of the data acquisition and processing applied regarding profiles SL1 and SL2, which were acquired in 2014, are reported in Polom et al. (2018). Profile SL6 was acquired in 2016 using the same acquisition method and parameters. The profile was positioned perpendicularly to the main structure dip of the Wadi Ibn Hamad alluvial fan with the intention of imaging a cross section of probable fluviatile subsurface channel structures at the eastern boundary of the alluvial fan.

3.2.4 Ground-penetrating radar

Ground-penetrating radar uses the active emission of electromagnetic waves in the microwave band from a transmitter to a receiver antenna to investigate the relative dielectric permittivity of the shallow subsurface (Fig. A1d). The physical background of GPR is explained in Appendix A.

All GPR surveys were conducted in October 2014 in common offset mode (transmitter and receiver antennas move simultaneously across a linear transect at a fixed distance) using a RAMAC GPR CU II system from Malå Geosciences (Fig. 4f). Several transects were collected using both 100 MHz and 50 MHz unshielded antennas, although only results from the latter are reported here for a total of two transects (GPR 1 and GPR 2 in Fig. 1b) that followed ERT surveys. Time-to-depth conversion for radargrams was based on estimates of electromagnetic wave velocities from three different approaches: (1) from diffractions in common offsets, yielding velocities ranging between 0.09 m ns⁻¹ for deeper materials and 0.15 m ns⁻¹ for surficial sediments, (2) from field calibrations using a metal rod at a depth of 1.2 m buried on the wall of a collapsed sinkhole, with values ranging between 0.15 and 0.2 m ns⁻¹, and (3) from samples collected in the field at approximately 2 m deep in a sinkhole wall and measured in the lab under different moisture contents, reaching average values of 0.15 m ns⁻¹ for the driest (5 % moisture content) conditions (Appendix A). The mean dielectric permittivity value hence lies around ε_r∼ 4, typical for dry sand.

Processing of common offset profiles was conducted with ReflexW v9.0 by Sandmeier Geophysical Research (Sandmeier, 2019). Steps were limited to (1) substract mean (dewow), (2) background removal, (3) manual gain, (4) bandpass frequency, (5) Kirchhof migration (based on an average velocity of 0.15 m ns⁻¹), and (6) topographic correction. In some instances (i.e. GPR2 transect) the presence of isolated point reflectors (Neal, 2004) resulted in the presence of diffraction hyperbolas in common offset profiles, allowing the construction of 2D models of velocity.

3.3 Inversion and numerical modelling techniques

3.3.1 Inversion of SP and ERT data

In this study, all SP anomalies were assumed to be generated by a 2D polarized sheet that is inclined and thin and that has an infinite extent perpendicular to the SP profile, in the following called “strips”. Therefore, the model parameters are the position of the centre of the strip along the profile (x_o), its depth (h), the inclination angle (α_i) and the half-width of the strip (a) and the polarization factor or electric dipole density (k). The particle swarm optimization (PSO) method was applied to determine values of those parameters (Monteiro Dos Santos, 2010). The PSO algorithm was first described by Kennedy and Eberhart (1995): it is a stochastic algorithm that simulates features of the social learning process as sharing information and evaluation of behaviours. The algorithm considers a community of N different models, and each model is updated taking into account its lowest rms achieved so far (called the “pbest model”) until the best model by any model of the community (called the “gbest model”) is obtained. The final solution will then be the gbest model achieved at the end of the iteration process.

The resistivity data were inverted by using the commercial program RES2DINV (GEOTOMO) v. 3.5 (Loke and Dahlin, 2002) that applies the Gauss–Newton inversion scheme. In a 2D model the Earth section below the acquired profile is divided into numerous rectangular cells. The objective of the inversion procedure is to vary the resistivity value of each cell in order to find a collective resistivity response that best matches the apparent resistivity measured. The forward problem of the Jacobian matrix value calculation is hereby addressed by a finite-element simulation, that is, to calculate the apparent resistivity response of a specific resistivity distribution. A non-linear L1 norm optimization method with smooth constraints is used to calculate the change in the resistivity of the model cells to minimize the difference between calculated and measured apparent resistivity (Loke et al., 2018; Sjödahl et al., 2008).

Figure 5 2D hydrogeological model and simulation time setup. The hydrogeological unit abbreviations are DS – Dead Sea, AS – alluvial sediments, LS – lacustrine sediments and CS – conduit system. Units are coloured by their horizontal hydraulic conductivity (kf_h). The text gives a description of the modelling approach between the initial (a) and final (b) states, the latter after 50 years of simulated DS regression. The simulated karst network evolution is presented as a zoomed-in view in the insets within each image. Top of alluvium is considered freely draining. The top centre of the model corresponds roughly to the location of the outflow point of canyon CM5; the profile follows the NW section of the AA^′ line of Fig. 3.

3.3.2 Hydrogeological modelling

A 2D hydrogeological forward model was developed by using Modflow v. 2005, Mt3DMS and Seawat v. 4 as integrated in the FloPy environment (Bakker et al., 2016; Lee, 2018) to understand the effects of the falling Dead Sea level and of the hypothesized development of karst-related conduits on the groundwater system. With this model of falling base level, we simulated the evolution of groundwater level and salinity distribution in the superficial aquifer system along a 4 km long transect perpendicular to the shoreline of the Dead Sea (i.e. along a roughly NW–SE-orientated profile following the NW section of profile AA^′ in Fig. 3, with the centre approximately at the head of canyon CM5). The model salinity distribution provided electrical resistivity values that were compared to those estimated from the ERT inversions. The 2D hydrogeological model setup and approach are presented in Fig. 5.

The lithological units simulated in the model are (1) alluvial fan sediments (“alluvium”), which comprise sand and gravels and which form a superficial aquifer, and (2) lacustrine sediments (“mud”), which comprise interbedded evaporites (such as gypsum, aragonite, calcite, and halite) and carbonates in clay to silt size (Khlaifat et al., 2010) and which are considered an aquiclude (Abelson et al., 2006; Siebert et al., 2014; Strey, 2014; Sachse, 2015). Geologically, these lithological units represent, respectively, the terrestrial and lacustrine facies of the Pleistocene Lisan and/or Holocene Ze'elim formations. Alternating layers of alluvium and lacustrine sediments have been simulated to take into account the local geology at Ghor Al-Haditha, where a silty-clay layer lying under the alluvium is reported in boreholes and has been imaged seismically at 40–80 m below the surface at a distance of roughly 1 km from the ∼ 1967 Dead Sea shoreline (El-Isa et al., 1995; Taqieddin et al., 2000; Polom et al., 2018). A similar geometric arrangement of aquifers and aquicludes is shown from boreholes for the western shore of the Dead Sea (cf. Yechieli et al., 2006).

The model dimensions are 4000 m by 400 m (Fig. 5), and the model is discretized into a grid of 1 m wide by 0.5 m high cells. The Dead Sea (with a density of 1240 g L⁻¹, e.g. Siebert et al., 2014) is included towards the west. Hydraulic heads are calculated based on a reference datum defined at 668 m below mean sea level. No flow is assigned to the bottom of the aquifer system. Boundary conditions of constant head are defined at the Dead Sea level in the north-west and a dropping water table by 0.75 m yr⁻¹ in the south-east, according to the nearest available well measurements (cf. Sect. 2). The shrinkage of the Dead Sea is simulated as a yearly fall in the hydrological base level by 0.5 m vertically for 10 years, 0.75 m vertically for 20 years and 1 m vertically for another 20 years according to estimations by Watson et al. (2019).

Table 2 Hydrogeological model material parameters.

The hydraulic properties in the model are mostly estimated in the absence of laboratory or in situ data, and they are chosen to represent a simple aquifer–aquiclude system (Table 2). The horizontal hydraulic conductivities of k_h = 1 × 10⁻⁶ m s⁻¹ for alluvium and k_h = 1 × 10⁻¹⁰ m s⁻¹ for lacustrine sediments are taken as mean values from the literature as adapted for the DS in studies from Shalev et al. (2006), Strey (2014), and Sachse (2015). An important assumption adapted to the local geologic findings is that vertical hydraulic conductivities are assumed to be 10 times lower than the horizontal ones. This represents anisotropy imparted from sedimentological layering and assumes the dominant role of bedding over tectonic-related anisotropy (i.e. fault-controlled flow) in this particular area (cf. Sect. 2). Hydraulic conductivities for the Dead Sea and for the hypothesized conduit system have been selected after various parameter tests. We added 1 m deep horizontal conduits of high effective porosity and hydraulic conductivity (Table 2) to the system in each fifth year of simulation. Due to the groundwater-level decline and shoreline recession, the horizontal extent of later-formed and deeper-lying conduits is greater than of earlier-formed and shallower conduits. Furthermore, we added drain cells to the whole surface of the alluvium cover to account for the surface channel systems that existed before regression of the Dead Sea (cf. Sect. 4.2.5). Conduit effective porosity and hydraulic conductivities were first varied in a trial-and-error approach to achieve model convergence and realistic steady-state hydraulic results, as further explained below in Sect. 4.3.

Table 3 Hydrogeological model simulation parameters.

We simulate density-driven groundwater flow based on salt concentration in the porewater. The modelling approach is to first simulate the saline water distribution in an estimated steady state over 100 years before the Dead Sea started retreating and to use this as the initial condition for simulating the yearly lake recession. The starting salt concentration of 340 g L⁻¹ is adopted for the Dead Sea and for the lacustrine sediments, which are assumed to be initially saturated with Dead Sea water. A constant salt concentration is defined for the Dead Sea only. “Salt” hereby stands as a proxy for all types of dissolved evaporite minerals in the numerical modelling, with NaCl as the main component. The lacustrine sediments are assumed to be subject to a continuous fresh or slightly brackish water inflow of initial salt concentration 0.66 g L⁻¹, through the alluvial sediments, which is based on values of the Ain Maghara spring (Table 1). A hydraulic head gradient of ∼ 30 m km⁻¹ (cf. Sect. 2; El-Isa et al., 1995; Sawarieh et al., 2000) is set, which is similar to the topographic gradient in this area (Al-Halbouni et al., 2017). Basic diffusion and advection processes are added by using appropriate parameters for the different aquifer–aquiclude materials (Table 3). Convergence limits of the changes in hydraulic head between two iterations were set to 1e-4; the same holds for changes in the concentration. Formation factors, effective porosities and the salt diffusion coefficient are known for the Dead Sea sediment analysed samples of the western side of the lake (Yechieli and Ronen, 1996; Ezersky and Frumkin, 2017). Transport steps and dispersivities have been derived from different saltwater intrusion studies and the classical Henry problem (Croucher and O'Sullivan, 1995; Geo Slope, 2004; Abarca et al., 2005; Langevin and Guo, 2006; Zidane et al., 2012; Meyer et al., 2019). For an overview of the initial spatial distribution of the effective porosity, concentration and head parameters, please refer to Appendix B.

Figure 6 Evolution of stream channels and karstic depressions from 2000 to 2017. Left column shows aerial or satellite imagery. The image sources are aerial survey of the RJGC (2000); Quickbird satellite (2006, 2012) and Pleiades satellite (2018). Right column shows maps of channels (red: straight/braided, green: meandering), ground cracks denoting the limits of a large-scale depression, and depression or sinkhole-hosted ponds (blue). Flow has converged upon CM5 in the years following 2012, with the other channels ceasing to prograde shoreward.

We estimate the bulk soil electric resistivities ρ in our models from the salt concentration by using transformation equations between chloride concentration and TEM resistivities (Ezersky and Frumkin, 2017). Hereby we rely on water resistivity vs. salinity relationships developed from borehole studies (Ezersky and Frumkin, 2017; Yechieli et al., 2001) and refer to the TDS content rather than the chloride concentration in our simulations.

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\text{TDS}=& {\mathit{\rho }}_{\mathrm{w}}^{-\mathrm{1.229}}\cdot \mathrm{3.019}\\ & \mathrm{\text{for low concentrations}}(\text{TDS}<\mathrm{131}\mathrm{g}{\mathrm{L}}^{-\mathrm{1}})\\ \text{TDS}=& {\mathit{\rho }}_{\mathrm{w}}^{-\mathrm{2.18}}\cdot \mathrm{0.247}\\ & \mathrm{\text{for high concentrations}}\left(\text{TDS}\mathrm{131}\text{–}\mathrm{320}\mathrm{g}{\mathrm{L}}^{-\mathrm{1}}\right)\end{array}\end{array}$$

We can derive the bulk soil resistivities via $F={\mathit{\rho }}_{\mathrm{w}}/\mathit{\rho }$ by using Archie's law (Archie, 1942) even for the clay-containing DS sediments, as they do not present a cohesion and cation exchange effect according to Ezersky and Frumkin (2017).

Figure 7 Evolution of the canyon–spring–sinkhole system close to the former mud factory at the head of CM5. (a, c, e) show aerial (balloon and drone, 2014–2016) or satellite imagery (2018). (b, d, f) show maps of the channel (yellow), surface water (blue) and sinkholes (dark orange). Filled stars indicate active springs; empty stars indicate previously active springs. The size of the filled stars indicates the approximate relative contribution of that spring to the flow downstream of the T junction between waters of the easterly and westerly springs.

4 Results

4.1 Remote-sensing and field surveys

The new observations presented here regarding the geomorphological evolution of Ghor Al-Haditha are focussed on the area around the site of the former Numeira mud factory (Fig. 6), which is coincident with the area in which the new geophysical surveys were conducted (see Fig. 1b for an overview). Here, stream channels of two distinct morphologies have been cut into the exposed evaporite-rich mud deposits of the former Dead Sea lake bed: meandering (CM) and straight (CS). The heads of all meandering channels have developed at spring points (in most cases, one per channel). Such springs lie either at the alluvium–mudflat boundary or within the mudflat deposits, and they originated at points lying over 100 m seaward of – and several metres below – the 1967 lake level. The heads of straight channels do not correspond to spring points; they initially developed some distance out on the mudflat, downslope from the termination of the active alluvial fan of the Wadi Mutayl.

As the shoreline has progressively retreated, both channel types have grown seaward, incising new channel segments into the lacustrine deposits of the former lake bed as the shoreline retreats over time. While the straight channels also show upstream growth (e.g. CS1–3 in Fig. 6), most meandering channels show little or no upstream growth (e.g. CM1–4 and CM6 in Fig. 6). Established sections of both channel types also widen progressively with time. From field observations, channel widening is commonly associated with fault-delimited slumping of the channel sides (Al-Halbouni et al., 2017). Both channel types deepen with time in response to the fall in base level; vertical incision is the primary fluvial erosive response to base-level fall. The lower sections of the straight channels are commonly braided and contain deposits ranging in size from sand to cobble. Deposit grain sizes within the meandering channels are mud to silt.

Sinkholes began developing in the area around the Numeira mud factory sometime after 1994, having first developed in the south of the study area in the mid-1980s (Sawarieh et al., 2000; Watson et al., 2019). By 2000, a cluster of water-filled holes lay along the alluvium–lacustrine boundary to the north-east of the factory site, while another cluster had formed on the alluvium about 750 m to the south of the factory (Fig. 6). Between 2008 and 2012, sinkhole development in this part of Ghor Al-Haditha accelerated dramatically and migrated from both the initial clusters toward the factory site. In close spatio-temporal association with sinkhole development was the formation of a gentle uvala-like depression (Fig. 6). This depression formed on a scale much larger than the individual sinkholes, and much of its perimeter is delimited by numerous ground cracks and faults (Watson et al., 2019). By 2012, the section of this uvala along the boundary between alluvium and lacustrine sediments was occupied by a small lake.

A remarkable meandering channel, whose development highlights spatio-temporal links between subsidence and groundwater outflow, is CM5. This formed in the summer of 2012 with its head initially located in the middle of the mudflat, about 750 m from the alluvium–mudflat boundary (Fig. 6). Headward channel incision progressed rapidly upstream over 3 months, in association with the drainage of the lake that had formed within the uvala (see Al-Halbouni et al., 2017, for details). Incision and headward erosion continued at a slower rate up to 2018. The main springs feeding CM5 all lie near the centre of the uvala, at or near its deepest points. Since the establishment of CM5 in 2012, the growth of nearby meandering channels CM1–4 has diminished markedly. Similar spring-fed meandering channels further south-west also seem to have ceased developing shortly before 2012 (Fig. 1b). Also, aside from one locality ∼ 500 m to the south-west, sinkhole development within the uvala since 2012 has been focussed within a 200 m radius of the head of CM5.

Figure 8 Field impressions of conduits and karstic groundwater springs at Ghor Al-Haditha. (a) Nested sinkhole in the alluvium (located in the southern part of the uvala) intersecting a further cavity or pipe. Photo: Damien Closson, 2008. (b) Sinkhole in lacustrine mudflats (located just inland of CM7; photo taken 2018) with pipe disappearing toward the top of the picture. (c) Conduit formed in lacustrine mudflat deposits near to the shoreward termination of CM1 (photo taken 2014). The flow has excavated a conduit in the weaker mud deposits below a more competent evaporite lamination. (d) Spring s₄ at CM5 at baseflow rate (estimated discharge: 2 × 10⁻⁴ m³ s⁻¹). The main stream is fed by several similar seeps within the depression at the channel head. Photo taken: 2018. (e) s₄ during high flow 3 d after an intense rainfall event (estimated discharge: 2 × 10⁻¹ m³ s⁻¹). Photo taken: 2015, modified after Al-Halbouni et al. (2017). (f) Close-up of s₄ during high flow. The nature of the sediment load varies from fine sediments to pebble-sized clasts (up to 5 cm in diameter). Photo taken: 2015, modified after Al-Halbouni et al. (2017).

More detailed evidence of the dynamic development of erosional and subsidence features around the head of CM5 during the period 2014–2018 is shown in Fig. 7. After its establishment in 2012 (Fig. 6), channel CM5 was fed by a complex suite of springs whose activity varied considerably between field campaigns in 2014–2018. In 2013, three springs, s₀, s₁, and s₂, fed the channel, with flow from s₂ rapidly eroding a new path to the main channel between June 2013 and September 2014. By the October 2015 field campaign, most of the flow into CM5 was now provided by a new spring to the south, s₄. Active collapses of the s₄ stream head were observed over the course of a few days, between 25 and 28 October 2015, when a large rainfall event occurred, associated with flash floods locally elsewhere. These stream head collapses occurred in conjunction with sinkhole openings along a line 20–30 m directly upslope of the stream head (Fig. 7). Flow rates near the channel head were around 2 × 10⁻¹ m³ s⁻¹, and coarse-sand to pebble-sized material was suspended in the flow (Fig. 8e–f; cf. video supplement). The other springs had either dried up or provided much reduced discharge to the channel. By December 2016, these sinkholes had been obliterated, as continued collapse and erosion produced a new section of the channel with concurrent landward migration of spring s₄, which continued to be the main source of flow within the channel. The nature of s₄ was by now a proliferation of seepage points in the channel head (Fig. 7). Discharge at one of these seeps was measured to be 2 × 10⁻⁴ m³ s⁻¹ in 2018, 3 orders of magnitude lower than during the high-flow event of 2015, and only fine sand and silt particles were suspended in the flow (Fig. 8d). This low discharge is also in line with qualitative observations of a low discharge rate at CM5 in 2014 and 2016 (as compared to the 2015 high-flow event). Between December 2016 and October 2018, new water-filled sinkholes of considerable diameter (∼ 30 m) had formed west of CM5. One of these holes was a new source of discharge to the canyon and is labelled s₅.

Figure 9 Summary of SP and ERT survey results at the head of CM5. (a) Map showing the arrangement of the SP electrodes in an array and the resulting heatmap of potentials (with reference to a fixed electrode at the south-eastern side of the array) overlain on an orthophoto from December 2016. The locations of the inverted profiles SP1–SP5 and their associated strip model anomalies are labelled, along with the modelled flow directions for the anomalies. Also marked are the locations of springs s₁, s₃, and s₄. (b) Electric resistivity distribution along ERT profile 3, which used a Wenner configuration with an electrode spacing of 1 m. Data inversion gave an rms error of 4.4 % after eight iterations. The observed surface material and the locations and flow directions of SP strip anomalies are labelled. (c) ERT profile 4, which used a Wenner configuration with an electrode spacing of 2 m. Data inversion gave an rms of 10 % after five iterations. The observed surface material and the locations and flow directions of SP strip anomalies are labelled. (d) SP lines 1, 2, 4 and 5 data and derived best inversion model results with an indication of rms error. (e) Graphical representation of polarized strip anomalies matching the modelled SP line data. The direction of water flow along a strip to produce the inverted anomalies is generally from minus (red) to plus (black). Due to noise in the data, the exact peaks in the modelled self-potential curves are not always reached. All results from inversion of SP3 have been omitted from all figures for clarity; to see these results, see Fig. A2 in Appendix A.

In situ EC measurements, as well as persistent localized vegetation growth around the ponds and in streams (Fig. 7), indicate the discharge of brackish groundwater from these springs. EC values from October 2015 for each of the springs s₁–s₄ range from 13 to 22 mS cm⁻¹. Repeated in situ measurements in November 2018 at s₂ and s₄ were slightly elevated as compared to 2015 but show similar values of 26 and 22 mS cm⁻¹. For context, note that we measured EC values of 100–220 mS cm⁻¹ in 2015 and 2018 at other springs (mostly sulfurous) emerging from the salt-dominated deposits further north within the Ghor Al-Haditha area (e.g. at the head of CM7 in Fig. 1b).

Although limited, there is some direct evidence of subsurface groundwater conduits. At the shoreward termination of CM1, subsurface conduits were found in 2014 (Fig. 8c) and have been noted for other channels to the north of the study area in Ghor Al-Haditha. These subsurface conduits tend to be associated with more competent evaporite horizons within the lacustrine sediments, which appear to form the roofs of the conduits, while the weaker clays are excavated beneath. Additionally, occasional “pipes” have been noticed in the base of sinkhole collapses in both the alluvial gravels and lacustrine deposits (Fig. 8a and b), implying that these materials are locally competent enough to sustain conduits, such as that inferred to have fed CM5 in 2015.

4.2 Geophysics

4.2.1 ERT and SP results at the head of active spring-fed channel CM5

As demonstrated in the previous section, the area around the head of stream channel CM5 has been a highly active zone of erosion and subsidence since its formation in 2012 (cf. Fig. 7; Al-Halbouni et al., 2017; Watson et al., 2019). These rapid geomorphological changes, occurring at the interface between alluvial cover and lacustrine sediments, made this area an obvious target for geophysical investigation of possible subsurface structures (i.e. cavities or conduits) in both materials. Consequently, near-surface geophysical surveys with the SP (array) and ERT (lines 3 and 4) methods were performed (Figs. 1 and 9).

Using the PSO 2dD inversion method (cf. Sect. 3.3.1), we extracted and inverted data for five profiles when traversing the nearest-neighbour interpolated SP array (SP1–5; cf. Fig. 9a and d). The model misfit rms ranges from 5 % to 15 %, while the SP data error has been estimated in the field to be 4.5 %. Due to the 2D nature of the inversions, there is also a ∼ 10 m horizontal error in the location of each anomaly (in x–y space; cf. Fig. 9a). From these inversions, we derived the dimensions and orientations of the inclined “strips” of varying depths, widths and inclinations which best represent the modelled curves. The location and extent of these strips and the inferred flow directions are illustrated in Fig. 9e. The length of each strip corresponds to the half-width of the matching anomaly in Fig. 9d. The direction of water flow along a strip to produce the inverted anomalies is generally from minus (red) to plus (black). It should be emphasized that the modelled strips are by no means the only possible arrangements able to explain each SP anomaly, particularly in regard to the depth of the modelled strips: a long, deep strip may produce the same anomaly as a short, narrow one. For clarity, the results of profile SP3 have been moved to Appendix A, which also contains the detailed geometric and electric parameters of each inversion anomaly (Table A1).

The SP survey produced a complex dataset, the main anomalies of which are two large positive (red) patches with a maximum potential of ∂V = 50 mV, where δV = V−V_ref (V_ref being the potential at the reference electrode). These positive anomalies correspond to the pair of upflow anomalies (v) and (vi) along SP2 within alluvial gravel deposits which have subsided and fractured considerably with the formation of the uvala there. The “strip” anomalies are proportional in magnitude to the size of each patch. Additional upflow anomalies (i) and (iii) are present along SP1 with corresponding positive anomalies of ∂V ∼ 10 mV in the array and a horizontal flow anomaly (ii). Strip (iii) results in a long-wavelength anomaly from −10 to 10 mV in the second half of SP1. It appears to be flowing toward the locations of the formerly active spring s₃. Other anomalies are recorded along SP4 and SP5, though the signals here are much weaker and noisier given the error margin, and hence the apparent “upslope” flow directions of some of these anomalies have to be treated carefully.

ERT lines 3 and 4 offer overlapping coverage with SP1 and SP5 and generally show decreasing resistivity with depth, stratified into two principal areas. The lower sections of both lines reveal an extremely low-resistivity layer (ρ_A = 1–5 Ωm) beneath more resistive zones. The top of this extremely low-resistivity layer lies at 8–10 m depth beneath the centre of ERT3, and in the case of ERT4 it forms a “wedge” in the north of the profile which is coincident at the surface with the lacustrine deposits of the former lake bed in the final section of the profile. The upper, more resistive zones of ERT3 and ERT4 show a range of resistivities from low values (ρ_A = 5–30 Ωm) indicated by light blue to light green colours, over middle values (ρ_A = 30–100 Ωm) indicated by dark green to yellow colours and high values (ρ_A = 100–500 Ωm) of reddish/violet colour. Patches of reddish/violet colours of very high resistivity (ρ_A > 500 Ωm) occur where the surface is highly fractured, with significant ground cracks. ERT4 shows a clear resistivity gradient from SSE to NNW also within the resistive zone. Although there is little apparent correlation between the SP flow anomalies and any significant resistivity anomalies in the ERT data, all the SP strips are modelled at depths above the very low-resistivity area within the upper, more resistive layers in a range of resistivity values between 15 and 100 Ωm.

Figure 10 ERT line 5 along the dry stream channel CM1 at the interface between alluvium and lacustrine sediments and associated surface manifestations. The survey used an inter-electrode distance of 5 m in a Wenner configuration. (a) Pleiades satellite image from April 2018 near canyon CM1 showing the profile location. (b) Electric resistivity distribution achieved by inversion of apparent resistivities. An rms of 7.5 % was achieved after six iterations. (c) Field photo from October 2018 looking NW along channel CM1, with contact between the alluvial fan deposits and the mudflat deposits of the former lake bed visible in the foreground near to the trees. Person for scale. The Lisan Peninsula and the Judean Mountains on the Dead Sea's western shore are visible in the background. (d) Hand specimen of the former lake-bed deposits from the bank of channel CM1; see part (c) for the location. Lens cap for scale. (e) Post-depositional carbonate concretions in the alluvial fan deposits; see part (a) for the location. Lens cap for scale. The minerals deposited are a mixture of aragonite, calcite, gypsum and clay minerals, all derived from the marls of the Lisan Formation. This degree of cementation is remarkable for the alluvium and is likely responsible for the elevated resistivity of the shallow subsurface at this location.

4.2.2 ERT results at inactive dry spring-fed channel CM1

To investigate the similarities and contrasts between presently and formerly active groundwater-fed stream channels, we performed an ERT traverse (line ERT5) which crosses the boundary between alluvium and lacustrine sediments next to channel CM1, which was dry in 2018 (although it was weakly active during field campaigns in the years 2014–2016). ERT5 reveals, similar to ERT lines 3 and 4, an extremely low-resistivity layer (ρ_A = 1–5 Ωm), which again coincides with surficial lacustrine mud deposits (Fig. 10). The ground above this layer in the south of the profile is more resistive, but the spatial distribution of these resistivities is not uniform, with a range of resistivities from low (ρ_A = 5–30 Ωm), indicated by light blue to light green colours, over middle (ρ_A = 30–100 Ωm) indicated by dark green to yellow colours and high values (ρ_A = 100–500 Ωm) by orange colours. Similar to ERT4, there is a clear resistivity gradient from south-east to north-west also within the resistive zone. Locally, very high-resistivity bodies (ρ_A > 500 Ωm) appearing as patches of reddish/violet colours in the profile can also be observed. The extremely low-resistivity layer is in horizontal contact with the more resistive zones in the centre (200 m) of the profile. At 150 m horizontal distance the depth to the low-resistivity layer is already 35 m. This zone of extremely low resistivity becomes more conductive towards the NW. Between 100 m and 140 m along the profile at depths below 10 m, there is a “tongue” of more resistive (green) material, which links gradually to the conductor beneath at around 30 m depth.

Field observations of the surface deposits around the ERT profile indicate that the differences in lithological properties between the lacustrine mud and the alluvial gravels are responsible for their contrasting responses to the injection of electrical current. The lacustrine mud exposed in the nearby stream channel walls (Fig. 10d) primarily comprises a light brown marl composed of carbonate and clay minerals, with some laminated horizons of evaporite deposits such as aragonite and gypsum (∼ 0.5 cm thick, white-coloured laminations). Idiomorphic halite crystals of up to 1 cm diameter also occur throughout this material. In contrast, the alluvial gravels and conglomerates are composed of mainly limestone and dolomite clasts with some basalt clasts (cf. El-Isa et al., 1995; Sawarieh et al., 2000), which are weakly to strongly cemented by minerals such as calcite, aragonite or gypsum. These minerals are dissolved and re-precipitated post-deposition as concretions on the clast surfaces and in the intervening pore spaces (cf. Bookman et al., 2004).

Figure 11 Comparison of shear wave reflection and ERT profiles. (a) Seismic line 1 (modified after Polom et al., 2018b). The inserted yellow trapezoidal area marks the shape of ERT1. (b) Seismic line 2 (modified after Polom et al., 2018b). The inserted yellow trapezoidal area marks the shape of ERT2. (c) Electric resistivity distribution along profile ERT1 achieved by an inversion of apparent resistivities in a Wenner–Schlumberger roll-along configuration with an rms error of 11.9 % after nine iterations. (d) Electric resistivity distribution along profile ERT2 achieved by an inversion of apparent resistivities in a Wenner–Schlumberger configuration with an rms error of 3.8 % after six iterations. Borehole lithologies of BH1 and BH2 (projected) are derived from El-Isa et al. (1995) and recoloured from Polom et al. (2018). (e) Orthophoto image from 2014 of the damaged section of road crossing the uvala along ERT1. Six main fractures, F1–F6, are identified. (f) Field photo from 2014 of the subsided old asphalt road surface, which is offset on fractures F1 and F5 and locally overlain by the aggregate. (g) Field photo from 2014 of the offset of the old asphalt road surface on fracture F6.

4.2.3 Comparison of shear wave reflection seismics and ERT results in the main sinkhole area

To provide further context to the previous geophysical studies performed further inland in the study area (El-Isa et al., 1995; Sawarieh et al., 2000; Polom et al., 2018) and to investigate possible changes to the subsurface hydrology using electrical methods, we performed two new ERT surveys on the alluvial fan, at the southern edge of the uvala (for locations, see Fig. 1). ERT lines 1 and 2 are coincident with seismic profile lines 1 and 2 as reported by Polom et al. (2018), respectively (Fig. 11a–d). The seismic profilings were carried out along asphalt-topped roads; the ERT surveys were carried out on the soil next to these roads. Profile 1 transects the margin of a main uvala-like depression in this area, while profile 2 crosses a subtler and narrower zone of subsidence that extends south-west from the uvala. Numerous sinkholes formed adjacent to these profiles between 2002 and 2010, mainly within the uvala and the narrow subsidence zone (Watson et al., 2019); those around profile 2 were filled in by local farmers. Throughout the surveying time period (2014–2018), subsidence and related ground cracks (small faults or fissures) remained apparent along both roads, but especially around the point where the road crosses the uvala at the north-eastern end of ERT1 (Fig. 11e–g).

In general, both ERT profiles are characterized by a vertical resistivity profile consisting of multiple distinct layers, which agrees well with previous studies. A distinctive lobe with an uneven basal surface, defined across most of ERT2 and the south-western–central part of ERT1, is characterized by high (ρ_A = 100–500 Ωm) resistivity values in the upper 30–40 m of the subsurface. This is underlain by a zone of low (ρ_A = 5–30 Ωm) to extremely low (ρ_A = 1–5 Ωm) resistivity values at depths below 40–50 m. The high-resistivity layer is overlain by a thin (5–10 m thick) layer of moderate (ρ_A = 30–100 Ωm) resistivity, corresponding to irrigated topsoil. This apparent stratification is in broad agreement with the borehole results of El-Isa et al. (1995): the middle high-resistivity layer corresponds to interbedded sands and coarser gravels, below lower resistivities correspond to sand and even lower resistivities represent the silt and clay layers detected only at the bases of the boreholes. The seismic data in these areas consist of reflectors dipping gently to the north-west, representing the top sets of the underlying alluvial fan system (Polom et al., 2018). Additionally, the resistivity values presented for ERT line 5 of Alrshdan (2012), ρ_A = 40–200 Ωm, are largely similar to those of the uppermost 25 m of ERT1. Similar to the results of ERT4 and ERT5 near the mudflat, there is a resistivity gradient from south-east (400 Ωm) to north-west (60 Ωm) visible in ERT2.

A clear link between subsurface, low-resistivity anomalies in the ERT data and surface subsidence features is imaged in the north-eastern section of ERT1, where the profile intersects the southern limits of the uvala. A striking area of reduced resistivities is apparent, with its centre around 120 m along-profile at around 20 m depth. These low resistivities (ρ_A < 10 Ωm) form a “blob” around 50 m wide and 20 m high, with a rapid gradient at its edge to a surrounding area of greatly elevated (ρ_A > 500 Ωm) resistivity. This region of elevated resistivity continues to the base of the profile here (Fig. 11c). The corresponding part of the seismic section consists of strongly scattered, chaotic reflection patterns. The low-resistivity anomaly is directly beneath the area of maximum subsidence along the road, as shown in the aerial and field photos in Fig. 11e–g. A smaller, vertically elongate low-resistivity anomaly 50 m along ERT1 is also visible. This anomaly may be corroborated by ERT data from line 4 of Alrshdan (2012), which overlaps the north-eastern end of ERT1 and shows a similar low-resistivity (ρ_A < < 10 Ωm) approximately 10 m diameter area at a similar depth (< 10 m).

Another significant anomaly in ERT1 is the strong horizontal resistivity gradient between the previously mentioned region of high (ρ_A > 500 Ωm) resistivity to the north-east and a lower (ρ_A = 10–20 Ωm) region to the south-west, around 250 m along-profile (Fig. 11c). This linear feature is also visible in the seismic section: an offset of ∼ 2 m between a strong reflector at ∼ 80 m depth (just below the vertical extent of ERT1), downthrowing to the north-east, is marked by blue arrows around 260 m along-profile. The VES survey results of El-Isa et al. (1995) in this area also found evidence of a vertical disturbance of this nature, which they attribute to the presence of a fault.

Figure 12 Evidence of water flow in sediments of the Wadi Ibn Hamad alluvial fan as deduced from aerial photos and reflection seismic profile 6 connected to stream channel CM5. (a) Coloured Corona satellite image from 1970 overlain on a 2018 Pleiades satellite image, depicting the former surface channel network in the alluvial fan deposits. Water content and vegetation appear in blue (modified after Al-Halbouni et al., 2017). Flow directions in the areas of elevated water content, as determined from the hydraulic gradient toward the Dead Sea, are highlighted by arrows to the right of the flow pathways. (b) Shear wave reflection seismic profile 6 along the Amman–Aqaba highway near Ghor Al-Haditha (cf. Fig. 2). (c) Structural interpretation of seismic profile 6; analysis is described in Polom et al. (2018). Structures associated with development of a former channel system are recognized in the central part. (d) Shear wave interval velocity overlain on reflection profile 6.

4.2.4 Shear wave reflection seismics along a section of the Amman–Aqaba highway

The hydrology of the Wadi Ibn Hamad alluvial fan has evolved significantly since the inception of base-level fall at the Dead Sea in the 1960s. This is highlighted in Fig. 12a, which shows the 1970 configuration of the wadi's delta system. The presence of vegetation and water at the surface (represented by blue colours on the map) highlights numerous distributary channels for water emanating from the wadi, two of which diverge from a fork a few hundred metres from the fan apex and run from there toward the former Numeira mud factory site (Al-Halbouni et al., 2017). The north-eastern end of ERT1 appears to cross one of these old channels, near the low-resistivity anomaly detected in the subsurface there (cf. Fig. 11c). After 1970, this surface water distributary system was significantly modified by the excavation of a single, straight exit channel for the Wadi Ibn Hamad and by the development of agricultural fields in a grid-like arrangement (Fig. 1b). An extensive area of dense vegetation directly at the former Dead Sea shoreline is nonetheless indicative of a sustained source of sub-surface freshwater around the Numeira mud factory site.

Seismic profile 6 runs along the Amman–Aqaba highway and crosses this former distributary channel close to its fork (Fig. 12a). The seismic section in Fig. 12b reveals a vertically extensive paleochannel system in the superficial aquifer as highlighted in Fig. 12c, with the first interpreted horizons occurring at around 20 m depth and the deepest at around 120 m depth. The derived velocities between 240 and 600 ms⁻¹ for these alluvial gravels (Fig. 12d) represent differing degrees of sediment compaction and/or cementation (cf. Polom et al., 2018). The axis of this interpreted palaeochannel system coincides with intersection of the 1970 distributary channel with the profile line. Synsedimentary faulting is visible to the north of the central axis of the paleochannel with downthrow to the south-west.

Figure 13 Water table estimation from GPR. (a) GPR profile 1 running NE–SW parallel to ERT1 reveals a horizontal reflector between 400 and 500 ns two-way travel time (∼ 30–34 m depth below the surface). Note the gap between 190 and 210 m which is due to a surface artificial concrete channel. (b) GPR profile 2 running NW–SE along line ERT2 revealing a similar reflector at around 500 ns two-way travel time (∼ 38 m depth below the surface). Note that GPR profiles have been topography corrected via linear interpolation only; i.e. the depression is not shown in GPR1. Characteristic reflections have been marked by arrows. (c) GPR profile 2 electromagnetic wave velocity distribution. A characteristic drop in wave velocity is seen at the interface between the saturated and unsaturated zones.

4.2.5 Water table inferred from GPR

GPR common offset profiles were processed according to Sect. 3.2.4 prioritizing the enhancement of deep reflectors in order to infer the location of the water table during data collection in October 2014 (Fig. 13a and b). Both profiles show the presence of shallow (0.3 m ns⁻¹) diffractions indicative of air reflections and multiples that may result from ground-coupling effects or ringing across the shallow-most layer (due to the higher conductivity when compared to the shallow high-resistivity zone below, Fig. 11c and d).

GPR profile 1 shows two major air reflections (i.e. hyperbolas, black arrows in Fig. 13a), which correspond to buildings near the end of the line and near the gap in data (surface artificial concrete channel located between 190 and 210 m along the profile). In most parts of the profile, an undulating near-horizontal reflector between 400 and 450 ns is imaged (blue arrows), which corresponds to approximately 30–34 m depth when using an average EM wave velocity of 0.15 m ns⁻¹ as inferred from fittings to subsurface hyperbolic diffractions and lab tests. This reflector is interpreted as the limit between the unsaturated and saturated zones. The undulating nature of this reflector matches with the resistivity interface in the corresponding parallel ERT1 survey (Fig. 11c), despite the latter being taken in 2018, and can be attributed to lateral contrasts in EM wave velocity.

GPR common offset profile 2 (Fig. 13b) corresponds to line ERT2 (cf. Fig. 1b) and also shows the presence of air reflections attributed to electrical posts placed parallel to the transect with consistent lateral spacing of 20 m intervals. Reflector facies for the first 8–10 m depth are again characterized by apparent multiples (i.e. ringing) that seem consistent with ERT line 2 results showing the top 10 m high electric conductivity layer. Other diffractions (between 10 and 25 m depth) indicate signal velocities around 0.14–0.16 m ns⁻¹, consistent with very dry material. This zone corresponds to the high-resistivity range found at a similar shallow depth in ERT2. There is the presence of a strong reflector (blue arrows) at around 500 ns two-way travel time (or ∼ 38 m depth) when using an average EM wave velocity of 0.15 m ns⁻¹. Fitting of hyperbolic diffractions in GPR2 are used to generate a 2D velocity model based on root-mean square velocities (Fig. 13c). The wave velocity shows a characteristic drop from 0.15–0.12 m ns⁻¹ to about 0.08–0.09 m ns⁻¹ at about 40–45 m depth that is consistent with a change in saturation and thus further supports interpretation of the deep reflector in Fig. 13b as the water table. The reflector is less undulating than in GPR1, which may indicate the absence of strong lateral variations. When compared to the first 280 m of ERT2, this depth corresponds to the transition zone between moderate and low resistivity (10–100 Ωm).

Furthermore, laboratory investigation of alluvial sample material with different gravimetric water content θ revealed a conductivity of 0.0027–0.45 S m⁻¹, i.e. a resistivity of ρ_A = 2.2–370 Ωm for saturated (θ ∼ 20 %) and dry, in situ, alluvium (θ∼ 5 %), respectively (Sect. A4, Appendix A).

4.3 Hydrogeological modelling

To study the effects of the falling Dead Sea level upon the subsurface hydrogeology in the study area, we applied the finite-difference groundwater flow and transport model including density-driven flow from the MODFLOW family, MODFLOW2005, MT3DMS and SEAWAT integrated within a FloPy environment (Bakker et al., 2016). We performed distributed simulation modelling of density-driven flow of the hypersaline Dead Sea water and fresh groundwater undergoing a base level drop of 1 m yr⁻¹ for 50 years after 100 years of hydrogeological steady-state conditions. The presence of subsurface conduits is simulated in the model by assuming an equivalent porous medium (EPM), i.e. assigning increased porosities to linear arrangements of cells, with new conduits added every 5 years. In our model space, realistic simulations (as compared to studies on the western shore and elsewhere) of the evolution of salt concentration and groundwater head under density-driven flow can only be achieved by including these simulated conduits in the model space.

Figure 14 Hydrogeological final model results incorporating a simulated conduit network. (a) Water head distribution. (b) Salt concentration with 150 and 30 g L⁻¹ contours, and (c) bulk soil electric resistivity. Left scale shows different Dead Sea regression phases with an indication of vertical drop. LS stands for lacustrine sediments; AS are the alluvium sediments; DS is the Dead Sea and CS stands for the conduit system. The grey areas correspond to the unsaturated zone – i.e. the area above the water table. See Fig. 15 for a more detailed view on flow, resisitivity and concentration.

The main criteria for evaluating the quality of the hydrogeological model solutions were (1) model convergence and (2) realistic hydraulic head development without singularities or jumps, which is expressed as a smooth water table as the model converges. Modelling the contact between alluvial and lacustrine sediments as an aquifer–aquiclude boundary with highly contrasting hydraulic conductivities did not yield model convergence. Only a drastic and unrealistic increase in the bulk hydraulic conductivity of the lacustrine sediments – to values equal to or greater than the alluvium – gave a realistic hydraulic head (see Appendix B). On the other hand, a local increase in the hydraulic conductivities in both alluvium and lacustrine sediments – i.e. the inclusion of the simulated conduit network – achieves both convergent solutions and realistic hydraulic heads. The final model results incorporating such a conduit network are shown in Fig. 14.

Head. With the regression of the Dead Sea, the hydraulic heads fall, and the water table declines (Fig. 14a), initially at the shoreline and then propagating further inland. This causes a strong gradient in the hydraulic head and water table that shifts further upslope during the regression. The groundwater level consequently decreases in the alluvial cover in a short time. Formerly saturated areas become dry, e.g. only after 20 years of regression (or ∼ 12.5 m of Dead Sea level fall) the alluvium surface becomes dry. At a location of x = 2500 m (i.e. c. 500 m from the former Dead Sea shoreline in 1960), the resulting water table drop during the simulated last 20 years of estimates to ∼ 24 m, and in the whole period of Dead Sea regression it amounts to ∼ 72 m.

Figure 15 Close-up views of the model alluvium–mud contact area. (a) Simulated bulk electric resistivity. (b) Flow direction, typical flow values and salt concentration. Contours of 150 and 30 g L⁻¹ are indicated. LS stands for lacustrine sediments, AS are the alluvial fan sediments and CS is the conduit system. Note the grid to assist interpretation.

Concentration. After the initiation period simulation of 100 years under steady-state conditions, the salinity of the system (Fig. 14b) is initially in a quasi-equilibrium with fresh and hypersaline water (ca. 50–350 g L⁻¹) domains separated by diffuse zone of brackish-saline water that is inclined landward. The freshwater–saline water system reacts relatively slowly to base-level fall because the diffusion of the saline water is a slow process. The adjustment of the interface is most pronounced in the shallow subsurface, above 468 $\mathrm{m}\mathrm{b}.\mathrm{s}.\mathrm{l}.$, where it is caused by the development of the simulated network of conduits. A rather salty water composition dominates the initial flow in the conduit system with concentrations between 50 and 150 g L⁻¹, with a subsequent freshening to brackish composition (ca. 10–50 g L⁻¹). The fresh–saline interface system builds a complex shape locally, and different levels of saline water can be encountered in a hypothetical vertical slice, but lateral variations also appear. A slight increase in salinity in the alluvial sediments at depth can be observed as a subtle eastward shift of the 30 g L⁻¹ salinity contour in the model below 518 $\mathrm{m}\mathrm{b}.\mathrm{s}.\mathrm{l}.$, and this is expressed also as a subtle decrease in resistivity below 468 m here (Fig. 15). Regionally the interface evolves into a landward-inclined shape of a Ghyben–Herzberg style interface.

Resistivity. The soil resistivity shows a duality of responses as the model develops, increasing in the near-surface lacustrine sediments due to inflow of fresher groundwater in the conduits, but decreasing gradually (Fig. 14c) at the diffusion front. In the shallow subsurface, a so-called conductive (or diffusion) front propagation towards inland by ∼ 100 m in 50 years can be observed in the model. The resistivity decrease is especially true along the alluvium–mud boundary, but also in the implemented channel water system ranging far landwards within the alluvial sediments. At deeper parts of the model, the resistivity adjusts only slightly to the small changes in concentration.

For comparison with resistivity patterns observed in our ERT surveys, Fig. 15a shows the simulated bulk soil resistivity results for saturated areas. In the shallow sub-surface of the model (above 468 $\mathrm{m}\mathrm{b}.\mathrm{s}.\mathrm{l}.$), the contact zone between lacustrine sediments and alluvium initially existed as a very conductive lower layer and very resistive upper layer. Upon initial diffusion, however, the resistivities are modulated to give a SE–NW gradient in resistivity that is especially pronounced in the alluvium. As the model evolves with base-level fall, the conduits appear usually as more conductive lines in most of the alluvium, especially further inland, and as less conductive lines in the lacustrine sediments. Important features of the numerical modelling that may assist in interpreting the ERT inversion results are therefore as follows.

1. Decrease in the alluvial fan resistivities from 65 Ωm to lower than 5 Ωm when subject to saline and highly saline water saturation, with a SE–NW gradient.

2. Very low to low resistivities of conduits within alluvium sediment of 0.5–10 Ωm.

3. Increase in the previously brine-saturated lacustrine sediment resistivities from 0.1 up to 1 Ωm due to refreshening, especially within conduits.

Flow. Initially, sub-horizontal flow of the order of ∼ 2 × 10⁻⁷ m s⁻¹ and slow diffusive deep flow of the order of ∼ 1 × 10⁻⁹ m s⁻¹ occur in the alluvial sediments (Fig. 15b). At the contact between the alluvial and lacustrine sediments (at a distance between 2000 and 3000 m from the shore) which also corresponds to the site of the uvala-like depression formed around CM5, a complex, multidirectional flow network is observed to form in the latter years. The flow is initially horizontal and evolves into a mixture of horizontal and vertical flow during the regression of the DS. Locally near the discharge into the Dead Sea, higher velocities up to ∼ 5 × 10⁻⁶ m s⁻¹ are noted. At the last stage, strong preferential flow of ∼ 3 × 10⁻⁵ m s⁻¹ dominates in the main conduits. Above the conduits, upward and downward flow patterns evolve of ∼ 4 × 10⁻⁷ m s⁻¹, below the conduit slower infiltration dominates.

5 Discussion

Here we first discuss the geophysical results and their interpretation and then the hydrogeological modelling. In view of these aspects and the insights into the evolution of surface geomorphology in the study area as derived from the remote-sensing data, we then propose a refined conceptual model for the hydrogeological evolution in the Ghor Al-Haditha sinkhole area.

5.1 Interpretation of geophysics in the context of hydrological evolution

In terms of interpreting ERT results, it is important to note that it is generally difficult to distinguish the electrical response derived from geological properties and that derived from porewater properties without well data. Apparent resistivity depends only on the electric conductivity (i.e. the number of free electrons). Different materials can have equal or similar conductivity, and so the material represented by a given apparent resistivity is always non-unique. Small variations in the portions of metal (e.g. Fe) or salt (e.g. Na, K), whether contained within the solid or liquid phases, can result in large variations in their electrical conductivity. Additionally, one has to bear in mind the principle of equivalence in electric and electromagnetic methods (Kirsch, 2006). An inserted stack of thin highly conductive layers within a low conductive background layer structure can lead to similar apparent resistivities in ERT to a thicker, less conductive layer. The conductance of both features would image similarly, and interpretation needs to be done carefully, on the basis of additional data, boreholes and conductivity measurements.

Consistent with former studies in the area (Frumkin et al., 2011; Alrshdan, 2012), our geophysical surveys on the alluvial fan (ERT1 and ERT2, Fig. 11) indicate a general decrease in resistivity with depth, broadly stratified into distinct regions. On the basis of logs of sinkhole walls nearby (Taqieddin et al., 2000; Sawarieh et al., 2000), the uppermost region at 0–10 m depth with middle range resistivities (30–100 Ωm), likely represents a well-irrigated topsoil overlying alluvium with a high silt and/or clay content. This grades downward into a highly resistive (100–500 Ωm) region that lies at 20–45 m depth. The basal surface of this layer has a “hummocky”, uneven geometry, particularly in ERT1 (Fig. 11c), which may represent initial sediment geometry, typical different saturation grades and/or salinization (Gonçalves et al., 2017; Farzamian et al., 2019a, b). This layer correlates with the sands and gravels logged in BH1 and BH2 by El-Isa et al. (1995). Below this, a less resistive (10–100 Ωm) layer extends to at least 70 m depth. VES data and models from 1994 (El-Isa et al., 1995), constrained partly by the borehole BH1, suggested a locally complex distribution of dry alluvium and alluvium with fresh, brackish, or saline water. If we consider the resistivity ranges ascribed to different groundwater conditions by El-Isa et al. (1995) to be correct, then we would interpret the highly resistive region at 10–40 m depth to be equivalent to dry alluvial fan deposits and the lower more conductive region to represent the presence of groundwater. At the base of ERT2, there is a small area of low (2–5 Ωm) resistivity, which may represent brackish groundwater within alluvium or salty mud.

Other studies at the DS coastal area have determined resistivity values for brine saturated materials between 0.25 and 0.4 Ωm for lacustrine sediments, and between 0.45 and 5.8 Ωm for alluvial fan sediments (Yechieli et al., 2001; Ezersky and Frumkin, 2017). The groundwater resistivity has been estimated by the same two studies in the Dead Sea coastal area to lie between 0.247 and 0.765 Ωm. From our numerical modelling, highly conductive parts (with resistivities lower than 5 Ωm) may relate to either brine to saline water-saturated alluvium and or brine to saline water-saturated lacustrine sediments, and hence a geological distinction based only on the ERT surveys alone is impossible. More specifically, the conductive feature in ERT line 1 below the uvala (Fig. 11c) may equally be a patch of water-saturated salty mud or salty alluvium; a distinction is not possible. However, the projected boreholes BH1 and BH2 helped to identify the major lithological units. This information and a clear transition between resistive and conductive areas in the profiles and in the models offer indications of the fresh–saline interface rather than lithological boundaries.

Overall, these ERT results for profiles 1 and 2 indicate an apparent decline of the level of the water table in the study area relative to the levels found by El-Isa et al. (1995) almost 24 years previously. GPR profiles 1 and 2 hereby help to constrain the limit between saturated and unsaturated ground of 30–38 m below the surface in October 2014. This limit, seen as a clear reflection in GPR, lies near the top of the less resistive layer estimated by ERT (10–100 Ωm) taken 4 years later, corresponding to the sand layer in the borehole logs of El-Isa et al. (1995). If we take the known water table depth (recorded in January 1995) in BH1 of 20.5 m (El-Isa et al., 1995, pp. 83–84) and estimate the water table depth in October 2018 from ERT to be 35–45 m, then a decrease in water table depth is obtained for this area of 15–25 m at a rate of 0.65–1.1 m yr⁻¹. This rate of water table fall is similar to the average rate of 0.75 m yr⁻¹ observed at the Mazra'a pumping station between 1960 and 2010 (Goode et al., 2013) and is in line with the rate of base-level fall. Declines of groundwater levels at rates of ∼ 0.25 and ∼ 0.9 m yr⁻¹ are reported from boreholes, respectively, at Tureibe (1985–2012) and Ein Gedi (2000–2014) on the western shore of the Dead Sea (Abelson et al., 2017).

The combined results from ERT and SP potentially reveal a complex water-flow system via conduits from the edge of the main uvala toward the vicinity of CM5, the main active spring at the border of the alluvium and lacustrine sediments. Shallow subsurface channels can be resolved and are indicated by lower resistivity than the surrounding, deformed alluvial fan deposits. One instance occurs at the edge of the uvala-like depression and below the damaged road on profile ERT1 (Fig. 11c; Sect. 4.2.4). Others are found in association with the SP anomalies upslope of the stream channel head, near the interpreted base of the alluvium and just above the underlying lacustrine deposits (Fig. 9b and c). The suggestion that these low-resistivity ERT anomalies represent flow in such a conduit network is consistent with EC measurements of the springs at the head of CM5 during the 2018 field campaign. These measurements, between 22 and 26 mS cm⁻¹, would equate to resistivities between 1 and 0.5 Ωm. Similar resistivities were measured both at the low-resistivity anomaly at the edge of the uvala, and at around 420 $\mathrm{m}\mathrm{b}.\mathrm{s}.\mathrm{l}.$ elevation along ERT3, which is the level of spring resurgence at CM5. A caveat is that ERT, even with a small electrode spacing as used in this study, cannot resolve the exact geometry of metre to sub-metre diameter flow pathways at depth.

The SP data indicate a dominant vertical upwards groundwater flow near the outflow point at CM5 but also horizontal (upslope) and downward directions. This is perhaps in part fracture-controlled locally, although the patterns of the SP anomalies resemble the shape and extension of the canyon head structures from aerial view. Also, in such a hydraulic setting of local artesian springs (cf. Al-Halbouni et al., 2017), an upslope part of the flow can develop locally at the alluvium–mud contact, as indicated by the SP results (Fig. 9) and by the simulated flow pattern (Fig. 15a). The dominantly upward flow features may indicate that piping plays a role in sinkhole formation in the vicinity of stream channel heads, such as the case at CM5 in 2015 (cf. video supplement).

5.2 Feasibility of the hydrogeological model

The nature of aquifers formed at a boundary between freshwater and saline water is critically dependent on the physical and chemical contrasts between saline water and freshwater. Fluid flow driven by density contrasts between saline water and freshwater is generally described by two mathematical approaches: (1) the “sharp-interface approximation” (Drabbe and Ghyben, 1889; Herzberg, 1901; Pool and Carrera, 2011), which assumes two fluids of constant density, no mixing and a pronounced saline–freshwater interface, and (2) the density-driven approach, which accounts for mixing of the two water types within a “transition zone” (Henry, 1964; Croucher and O'Sullivan, 1995; Simpson and Clement, 2004). Hydrogeological modelling of groundwater flow using both approaches has been previously undertaken for sites around the Dead Sea area (Yechieli, 2000; Salameh and El-Naser, 2000; Shalev et al., 2006; Strey, 2014; Sachse, 2015; Odeh et al., 2015; Alfaro et al., 2017) and provides sophisticated understanding of the flow systems and problems related to groundwater extraction. However, density-driven modelling with consideration of the salt concentration distribution has so far only been presented for the western DS coastal aquifer by Strey (2014), who focused on factors controlling model performance, by Shalev et al. (2006), who studied salt dissolution and sinkhole formation with vertical faults as preferred groundwater conduits, and by Levy et al. (2021), who studied effects of base level rise or fall on the fresh–saline interface.

The main contribution of the 2D hydrogeological modelling in this paper is a conceptual understanding of the hydraulic system at the Ghor Al-Haditha sinkhole area. The goal of the modelling is not the exact reproduction of this system, which is complex, heterogeneous, and subject to strong fluctuations of the groundwater table due to flash-flooding events, but rather to achieve converging results for a realistic head (water table) and salt concentration. As for any modelling, limitations exist that should be borne in mind when appraising it.

There are five main limitations to the modelling here. Firstly, we chose a simplified initial condition of salty lacustrine sediments and non-salty alluvium that is based on limited subsurface constraints. The Dead Sea is known for long-term and short-term fluctuations of its level, with different sedimentation, evaporation and erosion periods (Bartov, 2002; Neugebauer et al., 2015; Levy et al., 2020). The resulting depositional inter-fingering of alluvium and lacustrine sediments, a geometry known from boreholes on the western side of the Dead Sea, was thus chosen to approximate this. Secondly, we assume an initial state of a rather sharp and fixed fresh–saline boundary that only changes from the onset of regression. Due to the limitation of computer resources, diffusion processes that would have been present for a few 1000s of years while the sea level was rising slowly or fluctuating by a few metres (Bookman et al., 2004) have not been considered. Thirdly, there is a scarcity of long-term well data measurements for validation of the simulated heads. The SE model boundary condition has been adjusted to the nearest available well measurements, and the derived water levels are informed by the 1994 VES data interpretation (El-Isa et al., 1995; Sawarieh et al., 2000), the gradient of which has been included as a starting condition. Personal communications from local farmers confirm that some wells bearing water at approximately 15 m depth around 10 years ago and extending down to 30 m depth have since fallen dry. This magnitude of groundwater head decrease is predicted in the alluvial cover in our model. Fourthly, we did not simulate ephemeral streams, or temporary wadi discharge and groundwater movement, but rather adapted a continuous homogeneous inflow from the eastern boundary, which leads to freshwater accumulations at the model boundary between the conduits and the Dead Sea. At least partly this accumulation has been observed in nature as a freshwater lake that formed near CM5 between 2006 and 2012 (Al-Halbouni et al., 2017; Watson et al., 2019). Finally, although in homogeneous systems groundwater flow can be assumed perpendicular to the shoreline (effectively 2D), in complex 3D aquifers, groundwater flow directions are deviated by geological heterogeneities and hence angular (Meyer et al., 2018a, b). It is clear that we cannot capture irregular conduit pathways and turbulence with the used EPM approach for karst conduit modelling, the limitations of which have been widely documented (e.g. Kovács et al., 2005; Hartmann et al., 2014). However, it is not uncommon to use EPM in karst studies for general investigation and conceptual understanding of the flow system (cf. e.g. Scanlon et al., 2003; Abusaada and Sauter, 2013; Ghasemizadeh et al., 2015) and it is computationally less expensive in comparison to other techniques, e.g. a discrete fracture network; cf. Goldscheider and Drew (2007). Despite these limitations, therefore, the main results from hydrogeological modelling help to understand the observed geomorphological changes and to interpret the geophysical results.

A first main finding is that a conduit system (or other form of greatly increased hydraulic conductivity) needs to develop across the interface between alluvium and lacustrine sediments to provide a discharge of the enhanced water flow resulting from the higher head gradient. This conclusion is based on dozens of different tested models of whom only a selection is presented in Appendix B. The karst conduits develop at the boundary between initially relatively impermeable lacustrine sediments and relatively permeable alluvium to drain the system fast enough, and therefore should be present also in areas of similarly changing hydraulic conditions. In our model, we simulate a tiered configuration of conduits, with this geometry being tied to the sea level. Indeed, such a configuration is reported for some carbonate coastal karst aquifers under similar conditions of shoreline regression and falling regional hydrological base level (Fleury et al., 2007; Bakalowicz et al., 2008; Bakalowicz, 2015). At Ghor Al-Haditha, our surface observations indicate that such conduits can be generated by physical erosion (“piping”) of the weak alluvial and lacustrine materials, by chemical erosion of the evaporite component of the lacustrine deposits or by both processes. These processes are not explicitly simulated in the model, but incorporation of these and the related mechanical feedbacks could be the subject of future work (e.g. Romanov et al., 2020).

A second important contribution of the hydrogeological model is to aid the interpretation of the geophysical results. We used one set of transfer equations between TDS and electric conductivity specific to the Dead Sea materials (Yechieli et al., 2006; Ezersky and Frumkin, 2017). One must be careful therefore regarding overinterpretation of the resistivity values here. The approach is based on the definition of empirical parameters, the formation factors, a classical issue in hydrogeology. The influence of the material on the resistivity is clearly visible in the model, and this is due to the intrinsic dependency of the regional transfer formula (1) on formation factor and material type. However, the results are in the range of expectations from our ERT surveys, and therefore we consider our approach to be suitable, despite material heterogeneities.

The assumption of material mixture in the simulated conduits seems generally viable and leads to relatively low resistivities due to the low formation factor resulting from a high void space (Ali Garba et al., 2019). In nature, erosion would enlarge these conduits, such that they may become more electrically conductive or resistive than their surroundings, depending on whether they are infilled by air or water, on the salt content of the water, and on the nature of the surroundings. The lower resistivity of model conduits within alluvium (Fig. 15) supports our interpretation of low-resistivity anomalies in the ERT profiles (Figs. 9 and 11) as water-filled conduits.

Overall, our simulation of continuous salt diffusion predicts differing effects on bulk resistivity in differing parts of the hydrogeological system as base level falls. At shallow levels (Fig. 15), resistivity generally increases in the lacustrine sediments and uppermost alluvium near the alluvium–lacustrine boundary as groundwater salinity is reduced primarily due to freshwater advection via the conduits. This simulated process of aquifer freshening due to base-level fall is compatible with concurrent decreases in electrical conductivity, elevation of the fresh–saline interface, and groundwater head of ∼ 7 m observed on the western Dead Sea shore in the EG26 borehole at Ein Gedi from 2011 to 2018 (Levy et al., 2021). At deeper levels in the model, however, there is a subtle saline water intrusion and corresponding overall resistivity decrease within the alluvium due to salt diffusion. The model resistivity values at these deeper levels (< 15 Ωm) agree broadly with the low values given by ERT at depths > 50 m under profiles 1 and 2. These values might support an interpretation of invasion of alluvium by saline water there, or, because of interpretative ambiguity, they might support an interpretation of a layer of salty lacustrine sediments (cf. Polom et al., 2018). More broadly, this simulated process of saltwater intrusion may in principle explain an increase in electric conductivity concurrent with a groundwater-level fall of ∼ 12 m observed in the EG19 borehole at Ein Gedi from 2006 to 2012 (Abelson et al., 2017), although other processes (e.g. salt dissolution) may alternatively be responsible there.

5.3 Surface stream channels, subsidence, and relationship to subsurface conduits

Recent work at the Ze'elim fan sinkhole site on the western side of the Dead Sea has highlighted critical interactions between surface water and groundwater within the evaporite karst system (Avni et al., 2016; Shviro et al., 2017; Arav et al., 2020). This builds upon past qualitative or anecdotal links drawn between high-rainfall periods and sinkhole formation at other sites around the Dead Sea (e.g. Arkin and Gilat, 2000; Taqieddin et al., 2000). Significant rainfall events and related flash floods at Ze'elim have been shown to cause rapid sinkhole collapses and wider longer-term subsidence, with undersaturated floodwater being drained by newly formed sinkholes and then resurging, nearly saturated with salt, shoreward (downstream) of these holes. Surface water ingress thereby promotes the local development of sinkholes, which in turn feed surface water down into evaporite karst conduits, which undergo chemical erosion and destabilize and form more sinkholes in a self-accelerating manner. This interlinked development of stream channels, sinkholes, and evaporite karst conduits thus occurs in a manner akin to our understanding of classical limestone karst settings.

Our study at Ghor Al-Haditha offers a “long view” of similar processes that complements the abovementioned studies on the western shore of the Dead Sea. Our remote-sensing and geophysical data cannot resolve the effects of individual flood events as achieved at Ze'elim (Avni et al., 2016) but instead indicate a major long-term shift in the hydrogeological system that is reflected in the geomorphological development of the alluvial fan section of the Ghor Al-Haditha study area. Specifically, we observe a shift of spring-fed stream activity, as evidenced by channel incision, from the south-west to the north-east of the study area over the past 30 years with a focussing of groundwater discharge at CM5 from 2012 onward (Figs. 1b and 6). This is broadly coincident with a shift in the area of active subsidence from the south-west to the north-east also (Watson et al., 2019). These shifts have occurred despite artificial excavation and straightening of the Wadi Ibn Hamad in the 1980s to guide surface flash-flood water to the south-west. Such artificial drainage has enabled flash floods to bypass most of the main area of subsidence and sinkhole formation studied here (i.e. the uvala around the former factory site). Therefore, in contrast to observations at Ze'elim, it seems that surface water influx into the main subsidence area has been minimal. Only since 2006 has surface run-off from the Wadi Mutayl shown signs of being partly directed into the uvala (Fig. 6), though timing-wise this appears to be a consequence, rather than a cause, of uvala development. As at Ze'elim, the eventual ingress of surface water from the Wadi Mutayl into the uvala may have contributed to its accelerated development during 2006–2012 (cf. Avni et al., 2016). Nonetheless, we interpret the observed shifts in discharge location and subsidence around the main Ghor Al-Haditha alluvial fan as reflecting an overall re-routing of groundwater flow within the alluvial deposits and any potentially underlying lacustrine deposits.

Such subsurface flow re-routing is common in karst: active conduits regularly become constricted or blocked via sediment deposition or conduit collapse (e.g. Palmer, 1975; Despain and Stock, 2005; Šušteršič, 2006; Simms and Hunt, 2007; Plan et al., 2009; Farrant and Simms, 2011; Brook and Murphy, 2017). New conduits, often termed “floodwater diversion conduits” (Palmer, 1991), commonly develop to bypass such blockages. Another contributing factor in the north-eastward shift in groundwater resurgence with time at Ghor Al-Haditha may be the reaction of the hydraulic potential field to the uncovering of the bathymetry of the former lake bed, which is spatio-temporally variable across the study area. This is reflected in the rate of shoreline retreat, which is rapid in the southern part of the study area and is slowest in the vicinity of CM5, relative to the 1967 shoreline (cf. Fig. 3 of Watson et al., 2019).

Our study also complements recent work on the western shore by highlighting a similar short-term link between rainfall events, elevated stream discharge, conduit erosion and sinkhole collapses. During the high-flow event of 2015, discharge at the main spring (s₄) clearly emanated from a conduit within the lacustrine sediments, about 6 m below the alluvium–lacustrine interface (Fig. 8e and f; cf. video supplement). This altitude level is comparable to the level of suspected conduits as indicated by the SP and ERT results from 2018. In recent studies of similar processes on the Ze'elim fan (Avni et al., 2016; Arav et al., 2020), high-flow events were inferred to cause substantial erosion and rearrangement of the conduit network. This erosion seems to be largely chemical (i.e. related to subsurface salt dissolution), as evidenced by high TDS values and $\mathrm{Na}/\mathrm{Cl}$ ratios in the discharging springs downstream of the sinkholes there. In contrast, we observed relatively low electrical conductivities of the water at the resurging at CM5 in both high- and low-flow states (2015 and 2018, respectively), in line with historical measurements in the area. This observation and the high density of vegetation at the CM5 channel head suggest that the salt concentrations of the groundwater here are low, thus favouring mechanical erosion as a key agent of conduit development in this particular case.

These observations therefore highlight the potential role of physical erosion and the piping mechanism for generating conduits and sinkholes within the upper few metres of alluvial fan deposits and/or lacustrine deposits, as initially proposed by earlier studies at the Dead Sea (Arkin and Gilat, 2000; Taqieddin et al., 2000; Al-Halbouni et al., 2017). Further observations supporting the possibility of conduits being sustained in the alluvium include sustained water flow at the base of sinkholes in alluvium on the Neve Zohar fan on the western side (Arkin and Gilat, 2000). Such conduits are likely generated by the washing out of fines within the alluvium, especially during high baseflow following recharge of the fan by rainfall and/or floods, and they can be sustained by the cementation of the alluvium – especially the older deposits. It is stressed that such a physical mechanism of conduit generation likely goes hand in hand with conduit development by chemical erosion (i.e. salt dissolution) occurring at deeper levels, similar to what has been demonstrated at several sites on the western shore (Yechieli et al., 2006; Avni et al., 2016). In this view, chemical erosion and deformation at deeper levels could help to create the initial secondary porosity within the alluvium that is then further amplified by focussing of water flow in a feedback loop similar to that proposed for salt dissolution.

The spatial extent and dimensions of any conduit network at Ghor Al-Haditha remain speculative due to a lack of direct observation; however, some basic inferences can be derived regarding such a conduit network. Firstly, the subsurface development of any kind of branchwork geometry, as is often observed in limestone caves, is unlikely due to the absence of a network of interconnected subsurface fractures and of concentrated surface recharge points to the conduit network from surface depressions. As discussed by Watson et al. (2019), we regard dissolution by surface water as almost absent at Ghor Al-Haditha. In our study area, the intergranular nature of initial porosity, combined with the observed substantial discharge variations, suggests that the pattern of conduit development may be a spongework of connected, enlarged voids and pathways (Palmer, 1991, 2007; Klimchouck et al., 2000). No direct observations of spongework conduit zones have been made in evaporite karst areas, and they are uncommon in the case of limestone karst (Palmer, 1991, 2007). The locations where spongework cave passages have been observed are often areas of mixing of different water chemistries in rocks with high primary porosity; a well-documented example is the mixing of freshwater and seawater in coastal caves, which are often termed “flank margin caves”. Although flank margin caves are overall not very numerous, they have been mapped in karst areas across the globe, such as along the coastal carbonate platforms of the Bahamas (Mylroie and Carew, 1990), the coasts of the Yucatan Peninsula in Mexico (Back et al., 1984; Smart et al., 2006), and several coastal cave systems in Italy, Sicily and Sardinia (Ruggieri and De Waele, 2014; D'Angeli et al., 2015; Arriolabengoa et al., 2017). Further examples of corrosive mixing of different water chemistries to produce spongework caves occur in the Prichernomorsky Basin on the northern coast of the Black Sea (Klimchouk et al., 2012), where water derived from deep sources rich in sulfuric acid has mixed with shallower matrix flow, causing aggressive groundwater erosion under hypogenic conditions. We suspect that deep-seated groundwater sources rich in sulfur also contribute to the groundwater resurging at the alluvium–mudflat boundary at Ghor Al-Haditha: a number of the springs feeding the stream channels, including s₂ at CM5 (Fig. 7) but particularly those to the north of the studied area (such as that feeding CM7), have a highly sulfurous odour. It has been noted that other groundwater springs in the vicinity smelling of sulfur seem to derive from the Ram and Kurnub sandstone units (Khalil, 1992; IOH and APC, 1995), which would likewise imply a deep-seated source to these waters at Ghor Al-Haditha.

Figure 16 Conceptual model of the hydrogeological system at Ghor Al-Haditha. (a) Oblique aerial view sketch map (altitude of the eye is around 600 m) showing the geological units as outlined in Sect. 2, based on Khalil (1992), structural geology, catchments of the major wadis and areas of surface flow observed in channels and inferred subsurface flow in buried paleochannels (in the alluvial fan deposits) and in a conduit network (below the uvala). Note that the locations of sinkholes are schematic and not as mapped by Watson et al. (2019). Line X–Y is the location of the cross section in part (b), which shows the hydraulic structure of the subsurface as inferred from our geophysical studies and the hydrogeological modelling.

5.4 Conceptual model of hydrological and geomorphological processes

The culmination of our study is a refined conceptual model of the surface and subsurface hydrological processes which have developed at Ghor Al-Haditha since the onset of base-level fall in 1967 and how these processes link to the geomorphological evolution of the study area in space and time (Fig. 16). The fall in base level in the study area has led to the development of a strong hydraulic gradient and an input of fresher groundwater at decreasing elevations within the alluvial fan deposits. We propose that this recharge is derived primarily from flow within the Wadi Ibn Hamad and secondarily from the smaller Wadi Mutayl. Leakages from the bedrock aquifers may also play a role, especially in the northern part of the Ghor Al-Haditha sinkhole site (north of the area shown in Fig. 1b), where the prevalence of sulfurous springs suggests involvement of water from the Ram-Kurnub sandstone aquifer. We propose that groundwater flow within the alluvial fan is guided towards the Dead Sea by matrix flow initially. Buried channel systems may act as groundwater pathways due to the enhanced permeability afforded by the unconsolidated coarse sediments deposited within them. Exactly how deeply these paleochannels extend into the subsurface is uncertain and likely highly variable, though the one imaged by seismic line 6 extends to depths of at least 120 m below the present land surface (Fig. 12).

Upon encountering the lacustrine mud and evaporite deposits in the subsurface nearer the Dead Sea shoreline, there is a transition to more concentrated and complex groundwater flow within a network of conduits that is developed by both dissolution and physical erosion. Water flowing within the shallowest conduits is brackish in nature, as indicated by its enhanced electrical conductivity as imaged by ERT (Fig. 11c) and measured EC values at the springs feeding CM5 during the 2015 and 2018 field campaigns. The surveys presented in this paper indicate that such conduits exist between 5 and 30 m depth between the head of CM5 and the centre of the alluvial fan, though it is possible that the conduit network extends to greater depths than this. Erosion within this conduit network has produced a widespread area of surface subsidence which we term an uvala (see Watson et al., 2019, for a full discussion of the processes governing uvala and sinkhole formation at Ghor Al-Haditha). These conduits direct flow into surface stream channels via springs of spatially and temporally variable discharge. The stream channels grow during periods of activity of these springs and cease to grow when that section of the conduit network has become inactive. Spatio-temporally, the activity of surface streams fed by groundwater has shown a north-easterly shift over time (Figs. 1b and 6), which is likely to reflect changes within the conduit network.

Additionally, our modelling results indicate that a complex fresh–saline water interface builds up over different levels. The shallow subsurface is dominated by “advection” transport processes within the matrix and conduit network, while the dominant process occurring at depth is slower diffusion across a broad “interface zone” between saline Dead Sea water and the fresh and brackish water within the matrix and conduit network. These saltwater intrusions of resistivities ≤ 1 Ωm may extend remarkably more towards inland due to the high salt concentration of Dead Sea water, in line with magnetotelluric studies of the deep brine system from Meqbel et al. (2013).

To help elucidate the true nature of the subsurface hydrological system in the study area and its links to surface geomorphology, future work should include gathering of more comprehensive information on the subsurface lithology and continuous discharge data for the springs and channels. An understanding of the temporal characteristics of discharge at the springs feeding these surface channels is critical for developing an improved comprehension of the links between subsurface channelization and surface collapse at Ghor Al-Haditha. Further borehole drilling would likewise provide the opportunity to improve the ground truth of the geophysical surveys undertaken in the study area and to improve our understanding of the nature of erosion in the subsurface.

6 Summary and conclusions

This study supports the hypothesis that the subsurface hydrological flow regime at Ghor Al-Haditha contains elements of diffuse (matrix) flow and concentrated, heterogeneous conduit flow. Concentrated subsurface flow appears to occur across the study area at depths ranging between 5 and 30 m, bringing fresher groundwater to the constantly retreating Dead Sea shoreline.

At the former shoreline, the distribution of surface water flow is directed by subsurface conduit flow. During the studied period, surface stream channels (fed by springs occurring at the contact between alluvial gravels and lacustrine sediments of the former lake bed) dried up as subsurface flow re-routed to converge upon a new network of springs draining into a now-dominant central channel (CM5) which formed in 2012 at the centre of a major karstic depression hosting hundreds of sinkholes. ERT and SP data from the head of this channel resolve shallow subsurface anomalies inferred to represent a complex pattern of multidirectional water flow.

Further inland, new ERT surveys reveal features which we interpret as representing concentrated subsurface water flow at around 30 m depth directly beneath the outer limits of the major karstic depression. GPR and ERT data give constraints on the groundwater table which has fallen by 15–25 m within 23 years. Still further inland close to the major wadi which drains the mountains to the east, shear wave seismic reflection surveys reveal buried paleochannel deposits up to 120 m deep. We hypothesize that the paleochannel now acts as a pathway for groundwater flow within the alluvial superficial aquifer towards the distal parts of the Wadi Ibn Hamad fan delta in the locality of the now-dominant central channel CM5, where there are finer (and potentially soluble) materials. This enhanced groundwater flow ultimately promotes the formation of (and further flow concentration into) a network of conduits made by piping and/or dissolution.

Numerical simulations also suggest that conduit flow helps to maintain the deep interface between freshwater and saline water at the shoreline as the Dead Sea level falls over time. Realistic simulations of the evolution of salt concentration and groundwater head under density-driven flow can only be achieved by simulating an increasing number of conduits of high hydraulic conductivity in the model space as the Dead Sea level (and consequently the groundwater level) falls over time after starting the model run with diffuse flow alone. Furthermore, simulated electric resistivities and SE–NW resistivity gradients are in good agreement with results from ERT.

These combined lines of evidence shed light on the increasing influence of underground conduit flow across the study area in space and time as base-level fall has proceeded. This process leads to a complex flow pattern at hydrogeological boundaries and plays a key role in the formation and spatio-temporal distribution of canyons and surface depressions that have formed extensively across the study area, confirming the results of previous studies, and may also act to maintain density-driven flow as the hydrological base level of the Dead Sea falls over time.

Appendix A: Detailed background of the geophysical methods and additional results

We here present details of the different geophysical methods presented in Sec. 3.2 and additional results that support the research presented in the main text. Figure A1 gives a simplified overview of the survey methodology for each geophysical technique applied in this study.

Figure A1 Simplified sketches of the survey layouts used for the geophysical methods in this study. (a) Self-potential with fixed base configuration applied to an array. (b) Electric resistivity tomography using Wenner and Schlumberger electrode configurations. (c) Shear wave reflection seismic principle using a geophone (elastic wave receiver) spread. (d) Ground-penetrating radar survey of a 2D profile.

A1 Self-potential

The following background on the SP method is based mainly on comprehensive works from Vichabian and Morgan (2002), Jardani et al. (2007), Richards et al. (2010), Allègre et al. (2010), and Jouniaux and Ishido (2012). For more details, the reader is referred to these authors and references therein.

Although rarely used as a stand-alone method, SP (Fig. A1a) is the only geophysical method able to resolve groundwater flow in many circumstances and on a large scale. It may be used to determine the depth of the unsaturated zone under certain conditions and to locate drainage paths and flow tubes (Kaufmann and Romanov, 2016; Voytek et al., 2016). The main contributing physical–chemical mechanisms are of an electrochemical, thermoelectric, telluric, and electrokinetic nature. The desired streaming potential signal (without disturbing effects) results from groundwater movement, an electrokinetic coupling between groundwater and soil minerals, that creates the electric double layer (EDL; cf. e.g. Overbeek, 1952; Corwin, 1990; Revil et al., 1999a, b).

Clay minerals carry negative superficial electrical charges due to crystallographic defects. Attracted positive charges (cations) therefore stay at the surface of clay mineral aggregates. In the contact zone between groundwater and soil, a solution of such neutral assemblies of ions develops, but there exist different zones where the cations largely outnumber the anions attracted by the charge difference. Physically it can be described following the Stern–Gouy–Chapman (Revil and Linde, 2006) model which divides the EDL into a high cation concentration Stern layer, a diffuse double layer and a constant concentration layer. A zeta potential, or electrokinetic potential, is defined as the electrical potential between the solid mineral as a shear plane and a zero-potential surface (usually the liquid). This definition leads to a negative zeta potential for the movement of positive charges and forms the physical background for the calculation of the self-potential or streaming potential. We derive the governing main equations for streaming potential calculation in the following. The coupled flow equations for the saturated case are after Sill (1983).

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}& J_{\mathrm{e}}=-{\mathit{\sigma }}_{\mathrm{0}}\mathrm{\nabla }V-L_{\text{ek}}\mathrm{\nabla }P\\ & J_{\mathrm{f}}=-L_{\text{ek}}\mathrm{\nabla }V-k_{\mathrm{0}}/{\mathit{\eta }}_{\mathrm{f}}\mathrm{\nabla }P\end{array}\end{array}$$

J_f is hereby the primary flow of matter (e.g. Darcy flow of water), and the effects of the coupled secondary flow J_e are usually considered to be negligible. However, here J_e is the point of interest. V is the streaming potential for constant temperature and no concentration gradients measured as a difference between the potentials in the liquid and in the grain. L_ek is the electrokinetic coupling coefficient following Onsager's reciprocal relation in steady-state conditions. η_f is the dynamic viscosity of the fluid, σ₀ the fluids' specific conductivity, k₀ the hydraulic conductivity and ∇P the pressure gradient.

Conservation of current (∇⋅J_e = 0), the electric convection current being compensated by the conduction current, leads to Poisson's equation with a source term:

$$\begin{array}{}\text{(A2)}& \mathrm{\nabla }\cdot \left({\mathit{\sigma }}_{\mathrm{0}}\mathrm{\nabla }V\right)=\mathrm{\nabla }\cdot \left(L_{\text{ek}}\mathrm{\nabla }P\right).\end{array}$$

For a homogeneous medium,

$$\begin{array}{}\text{(A3)}& {\mathrm{\nabla }}^{\mathrm{2}}V=C{\mathrm{\nabla }}^{\mathrm{2}}P,\end{array}$$

with C as the streaming potential coupling coefficient defined when J_e=0:

$$\begin{array}{}\text{(A4)}& C=-L_{\text{ek}}/{\mathit{\sigma }}_{\mathrm{0}}=\partial V/\partial P.\end{array}$$

The streaming potential coupling coefficient for unidirectional flow of a liquid (anode) between two electrodes (cathodes) derived from the classical electroosmotic model (Overbeek, 1952) can be expressed as

$$\begin{array}{}\text{(A5)}& C=\partial V/\partial P={\mathit{\epsilon }}_{\mathrm{f}}\mathit{\zeta }/{\mathit{\eta }}_{\mathrm{f}}{\mathit{\sigma }}_{\text{eff}},\end{array}$$

with σ_eff=Fσ₀ as the effective conductivity and ε_f the permittivity of the fluid. ζ is the zeta potential inside the electrical double layer. This Helmholtz–Smoluchowksi equation is valid for a small double-layer thickness (Debye length) in comparison with the pores and capillaries and at constant viscosity and permittivity if fluid conductivity dominates. Typical SP gradients of groundwater flux along an axis x are around $\mathrm{d}V/\mathrm{d}x$ = 0.2 mV m⁻¹. For karst aquifers this equation holds for any capillary and the SP signal comes from capillaries and the matrix. A highly amplified signal is created if e.g. sand is in the fracture. The signal is a combination of topographic regional-scale flow and local heterogeneities imposed by the karst, e.g. sinkhole downward seepage.

Figure A2 SP data and inversion results for profile SP3.

The detailed results of the inversion of five SP profiles extracted from the SP array are summarized in the following table, and the result for SP3 is presented in Fig. A2. They represent the model parameters of geometry and electric charge of the modelled thin sheets (called strips) that extend infinitely horizontally perpendicular to each profile. The shape factor q approaches zero for such modelled horizontal strip-like structures. The dipole moment for inversion has been constrained to positive values only and relates the current density J to resistivity for a thin polarized strip-like structure to $K=J\mathit{\pi }{\mathit{\rho }}_{\mathrm{A}}/\mathrm{2}$. For more details on SP modelling and inversion, see Murthy and Haricharan (1985), Monteiro Dos Santos (2010), and Biswas (2017).

Table A1 Inversion results for SP modelling with the PSO method of Monteiro Dos Santos (2010).

A2 Resistivity tomography

ERT works by injecting direct-current (DC) or low-frequency alternating current (AC, short-time DC application) into the ground and measuring the differences of the electric potential field at multiple electrodes (Fig. A1b) along a profile or array via computerized tomography (see e.g. Kirsch, 2006). The resulting apparent resistivity is based on Ohm's law, and the true resistivity of the ground is calculated thoroughly by 2D or 3D inversion (e.g. Dahlin and Loke, 1998; Loke and Dahlin, 2002; Jardani et al., 2013; Loke et al., 2018). Small electrode spacing leads to high resolution and shallow penetration depths, whereas large electrode spacing lowers resolution and increases penetration depths.

Measuring the apparent resistivity to ρ_A = $K^{U/I}$ by electric current I and voltage U needs the estimation of geometric factors K for each survey layout. Generally, these are calculated by the distance differences between two electrodes (A, B, M, N):

$$\begin{array}{}\text{(A6)}& K=\mathrm{2}\mathit{\pi }{\left[\left({\displaystyle \frac{\mathrm{1}}{AM}}-BM\right)-\left({\displaystyle \frac{\mathrm{1}}{AN}}-{\displaystyle \frac{\mathrm{1}}{BN}}\right)\right]}^{-\mathrm{1}}.\end{array}$$

For the Schlumberger configuration we get

$$\begin{array}{}\text{(A7)}& K=\mathit{\pi }\left({\displaystyle \frac{a^{\mathrm{2}}}{b}}-{\displaystyle \frac{b}{\mathrm{4}}}\right).\end{array}$$

For the Wenner Alpha configuration we get

$$\begin{array}{}\text{(A8)}& K=\mathrm{2}\mathit{\pi }a.\end{array}$$

Estimation of the penetration depth however is difficult subject to the concept of conductance vs. thickness of a layer, as both can influence the distribution in the same way. The integrated conductivity can be derived by inversion methods as described by Loke et al. (2018) and in the corresponding part of the paper (Sect. 3.3.1).

A3 Shear wave reflection seismics

A seismic source at or close to the Earth's surface generates an elastic wave signal that propagates into the subsurface. The seismic wave is partly reflected at interfaces between elastically contrasting layers and then recorded by receivers (geophones) at the surface (Fig. A1c). Another part of the wave energy is transmitted through the first layer interface and reflected and transmitted at the next layer interface. This continuous process enables the imaging of stacked subsurface-layer structures with superior resolution and depth penetration compared to other geophysical methods. The main factors influencing the structure imaging are the wave propagation geometry, the propagation velocity and the signal propagation time from the source to the receiver. Shear body waves (S waves) are not affected by ore space content (fluids, gases), and they propagate significantly more slowly than compressional body waves (P waves). The resulting smaller wavelengths of S waves can improve the subsurface resolution by up to 10 times. Furthermore, the S waves propagate without any disturbance of groundwater, which is advantageous in alluvial sedimentary systems like the Wadi Ibn Hamad delta fan.

A4 Ground-penetrating radar and sample laboratory tests

The GPR method typically achieves vertical resolutions between 1 and 10 s of centimetres depending on antenna frequency (Jol, 2009). The method has been widely used for decades in various sedimentological characterization studies (e.g. Neal, 2004) and has proven successful in identifying sinkholes and when applied to more general characterization of karst systems, often in combination with other geophysical methods (e.g. Kruse et al., 2006; Gutiérrez et al., 2011; Gómez-Ortiz and Martín-Crespo, 2012; Carbonel et al., 2014; Mount et al., 2014; Margiotta et al., 2016; Sevil et al., 2017; Kaufmann et al., 2018; Čeru and Gosar, 2019). The method has proven especially useful when determining the limits between unsaturated and saturated groundwater zones, imaging the water table in karst systems, and defining the freshwater–saltwater boundaries in coastal aquifers (e.g. Kruse et al., 2000; Gustafsson, 2005; Margiotta et al., 2012; Mount et al., 2015; Paz et al., 2017).

An electromagnetic impulse of a distinct frequency f is sent via a transmitter antenna into the ground, and the ground response waves are recorded by a receiver antenna mounted usually close to the sender (the so-called vertical incidence case, Figure A1d). The survey is performed by dragging the ground-coupled antennas along the surface with a steady movement to generate regular measurement intervals and setting digital marker points at distinct horizontal distances to later match resulting trace numbers (individual measurements) against horizontal distance.

Figure A3 Results of electromagnetic laboratory tests on alluvium samples for three different moistures θ. (a) Real part of dielectric permittivity ε. (b) Imaginary part of dielectric permittivity ε. (c) Direct-current conductivity σ_dc. (d) Attenuation α. The results show that attenuation in the moist alluvial material is high due to high DC conductivity. For dry alluvium damping is comparably low.

In GPR surveys, water content may lead to strong damping of the electromagnetic signals by direct current losses, preventing sufficient depth penetration. Therefore, laboratory measurements have been performed to determine the properties of alluvial material samples from the field site. Figure A3 contains the results for electromagnetic laboratory tests on the samples with a frequency of f = 100 MHz. Samples had a high iron mineral and salt content, and, together with gravimetric soil moisture θ between 5 (in situ) and 20 %, representative for this irrigated area, the conductivity σ_dc = $\mathrm{1}/{\mathit{\rho }}_{\mathrm{A}}$ lies between 2.7 and 0.5 S m⁻¹. The relative dielectric permittivity ε lies between 4 and 24 and also depends on the moisture content. This results in an attenuation of 10–100 dB m⁻¹ for 100 MHz, and the penetration depth δ, or skin depth, for a 100 MHz antenna in a homogeneous half space can be approximated to δ = $\sqrt{\frac{\mathrm{2}{\mathit{\rho }}_{\mathrm{A}}}{\mathit{\omega }\mathit{\mu }}}$ (Jol, 2009) with μ = μ_rμ₀ = $\sim \mathrm{1}\cdot \mathrm{12.57}\times {\mathrm{10}}^{-\mathrm{7}}$ N A⁻² as the magnetic permeability and ω=2πf as the angular frequency. This results in δ∼ 35 cm for a wet soil (σ_dc = 0.45 S m⁻¹). For dry soil (σ_dc = 0.0027 S m⁻¹), the penetration depth would be high (∼ 970 m) if no other attenuation effects were apparent.

The signal frequency also has an impact on penetration depth and resolution of the GPR surveys (Jol, 2009). A low-frequency antenna was preferred in this survey because the target, void structures and groundwater table are expected to be at more than 10 m depth, although this incurs a low horizontal resolution of Δl = $\sqrt{\frac{v_{\text{em}}r}{\mathrm{2}{f}_{\mathrm{c}}}}$ with r as the radial distance from the sender. For 50 MHz and 0.15 m ns⁻¹ velocity, the minimum horizontal extension of an object that can be resolved at 10 m depth is 3.87 m. For 100 MHz it is 2.7 m.

Appendix B: Hydrogeological model parameters, alternative models and resistivity calculation

The hydrogeological modelling procedure as described in Sect. 3.3 requires the definition of further important parameters such as effective porosity and concentration distribution, initial conditions and fixed head and concentration locations (Fig. B1).

Figure B1 Hydrogeological model parameters. (a) Initial effective porosities for all materials. (b) Location of fixed head (−1), otherwise free head (1). Note that boundary conditions are also fixed but not visible due to the small cell width. (c) Initial location of fixed concentration (−1), otherwise free concentration (1). At each channel mouth the DS concentration is set to free adjustment. (d) Initial salt concentration distribution. Dead Sea and lacustrine materials are valued by high values (340 g L⁻¹), the alluvial sediments by low values (0.66 g L⁻¹). LS stands for lacustrine sediments; AS are the alluvium sediments; DS is the Dead Sea and CS stands for the conduit system.

Alternative models for a convergent steady-state solution have been tried and are presented in Fig. B2. There are no channels installed in these models; only the hydraulic conductivities of the lacustrine sediments are varied. The water table does not follow a realistic curve in either alternative model 1 or alternative model 2. In alternative model 1 we achieve a more realistic shape of the water table, but the overall magnitude of the hydraulic conductivity of LM had to be increased to at least the same order of magnitude as the AS material, an unrealistic scenario for the whole lacustrine material. Hence, these models were discarded in favour of the model presented in Sect. 4.3 with local conduits of high permeability.

Figure B2 Alternative (deprecated) hydrogeological model results for concentration and water table. (a) Model 1 with horizontal and vertical hydraulic conductivity of LS of the same value as AS (k_h = 1 × 10⁻⁵ m s⁻¹, k_v = 1 × 10⁻⁶ m s⁻¹). This high hydraulic conductivity for mud is unrealistic for the initial aquiclude conditions. (b) Model 2 with horizontal and vertical hydraulic conductivity of LS 4 orders of magnitude lower (k_h = 1 × 10⁻⁹ m s⁻¹, k_v = 1 × 10⁻¹⁰ m s⁻¹) than of AS, like in the final model but with no karst conduit system. In both examples the water table behaves in an unrealistic (non-linear) way, and non-convergence and no clear saltwater–freshwater boundary changes are observed.

To calculate the resistivity from salinity simulations, one can alternatively to Sect. 3.3.2 use the direct bulk soil resistivity–salinity relationships from Ezersky and Frumkin (2017).

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}& \text{TDS}={\mathit{\rho }}_{\mathrm{A}}^{-\mathrm{3.35}}\cdot \mathrm{15.5}\\ & \mathrm{\text{for high concentrations}}(\text{TDS}\ge \mathrm{131}\mathrm{g}{\mathrm{L}}^{-\mathrm{1}})\text{in AS}\\ & \text{TDS}={\mathit{\rho }}_{\mathrm{A}}^{-\mathrm{1.377}}\cdot \mathrm{42.66}\\ & \mathrm{\text{for low concentrations}}(\text{TDS}<\mathrm{131}\mathrm{g}{\mathrm{L}}^{-\mathrm{1}})\text{in AS}\\ & \text{TDS}={\mathit{\rho }}_{\mathrm{A}}^{-\mathrm{2.61}}\cdot \mathrm{8.52}\text{in LM}\end{array}\end{array}$$

The results differ only in the magnitude by a constant shift to lower values of roughly 30 Ωm and are not shown here.