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

Karst groundwater systems are characterised 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 25 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 30 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 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 to represent subsurface water flow. Low-frequency GPR surveys reveal the limit 35 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 under-saturated groundwaters in response to base level fall perform realistically if they include the generation of karst conduits near the shoreline. The 40 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. https://doi.org/10.5194/hess-2021-37 Preprint. Discussion started: 1 February 2021 c © Author(s) 2021. CC BY 4.0 License.


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 (Goldscheider, 2015;Hartmann et al., 2014;Price, 2013). Consequently, a characteristically well-evolved network 5 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 (Hartmann et al., 2014;Kaufmann and Braun, 2000).
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 10 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 (Sauro, 2012;Waltham et al., 2005). 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 (Ford and Williams, 2007;Palmer, 1991Palmer, , 15 2007Palmer, , 2012. 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 20 semi-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). Base flow in these conduits may be of very low discharge, and surface flow may only occur rarely in ephemeral wadis (dry river valleys), where 25 water often flows beneath the surface (Price, 2013;Salameh et al., 2018;US Geological Survey et al., 2013).
One of the most rapidly developing evaporite karst systems on the globe is found on the shores of the Dead Sea ( Figure 1a). This hypersaline terminal lake, fed primarily by the River Jordan, 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 30 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 m msl to -434 m msl (1967ISRAMAR, 2020), i.e. by 40 m as of 2020. The lake level fell by 0.5 m yr -1 in the 1970's and by 1.1 m yr -1 in the last decade. This has led to a dynamic reaction of the hydrogeological system and landscape 35 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 lakebed that have been revealed by the retreat of the shoreline (see e.g. Abelson et al., 2006Abelson et al., , 2017Al-Halbouni et al., 2017;Watson et al., 2019;Yechieli https://doi.org/10.5194/hess-2021-37 Preprint. Discussion started: 1 February 2021 c Author(s) 2021. CC BY 4.0 License. et al., 2016 and references therein). The hazards posed by erosion and subsidence results in severe consequences for tourism, agriculture and infrastructure at the Dead Sea (Abou Karaki et al., 2019;Arkin and Gilat, 2000;Fiaschi et al., 2018).
There are several hypotheses regarding the nature of the subsurface hydrogeology at the Dead Sea and how this 5 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 increase promotes increased 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 lakebed (Bowman et al., 2007(Bowman et al., , 2010; Dente et al., 2017Dente et al., , 2019Ben Moshe et al., 2008;Vachtman and Laronne, 2013). Several studies contend 10 also that the base-level fall results in lateral migration of the fresh-saltwater interfaces in the subsurface, which enables relatively understaturated 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-15 conductive regional tectonic faults (Abelson et al., 2003;Charrach, 2018;Closson et al., 2005;Shalev et al., 2006). 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 (Kafri et al., 2007;Salameh and El-Naser, 2000;Yechieli, 2000) and (2) a density-driven flow simulations that account for mixing of waters of contrasting salinities and densities (e.g. 20 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 25 (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 (Abelson et al., 2017;Ezersky and Frumkin, 2013;Shalev et al., 2006;Yechieli et al., 2016). Other studies emphasise 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 ( Arkin 30 and Gilat, 2000;Sawarieh et al., 2000;Taqieddin et al., 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 generates subsurface cavities and conduits (Al-Halbouni et al., 2017;Arkin and Gilat, 2000;Polom et al., 2018;). It has been shown that under the geologic setting of the Eastern Dead Sea these void spaces are mechanically instable and may lead to different scales of 35 morphological expressions of subsidence . It is emphasised here that such varying views are not mutually exclusive (cf. Watson et al., 2019). Indeed some recent studies have emphasised linkages between the surface water flow, groundwater flow in a shallow conduit system and dissolution of deeper-https://doi.org /10.5194/hess-2021-37 Preprint. Discussion started: 1 February 2021 c Author(s) 2021. CC BY 4.0 License. seated salt -particuary during flash flood events -as an important factor in sinkhole formation (Arav et al., 2020;Avni et al., 2016).
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 5 (Goldscheider, 2015;Goldscheider and Drew, 2007). An overview of indirect methods applied to the Dead Sea sinkhole phenomenon is given by  and Salem (2020). Despite these challenges, near-surface geophysical methods have shown potential for cost-effective investigation of subsurface erosion in evaporite karst environments (see e.g. Abelson et al., 2018;Fabregat et al., 2017;Giampaolo et al., 2016;Gutiérrez et al., 2011;Jardani et al., 2007;Krawczyk et al., 2012;Malehmir et al., 2016;Muzirafuti et al., 2020;Wust-Bloch and Joswig, 10 2006) and determining the actual configuration of the fresh-saline water interface (e.g. Kafri et al., 2008;Kafri and Goldman, 2005;Kruse et al., 2000;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 for dissolution of a massive salt layer in the shallow subsurface (as described above) with cavity 15 formation as imaged by high geoelectric resistivities, and for linkage of this processes to sinkhole development in many areas there (cf. Ezersky, 2008;Ezersky et al., 2011;Yechieli et al., 2006). 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 20 on the Dead Sea's eastern shoreline, at a key part of the Ghor Al-Haditha sinkhole site in Jordan (Figure 1). We combine surface manifestations, remote sensing analysis, and multiple geophysical surveying methods of horizontally polarized shear wave (SH) 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. 25 First, we summarise 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 30 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 complements 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 hypothesise that flow of subsurface groundwater undergoes a transition 35 from 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 conditions in the case of a regressing Dead Sea level. This multi-disciplinary analysis https://doi.org /10.5194/hess-2021-37 Preprint. Discussion started: 1 February 2021 c Author(s) 2021. CC BY 4.0 License. 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.

Site of investigation and previous geophysical and hydrogeological studies there 15
The Ghor Al-Haditha study area is situated on the eastern shore of the Dead Sea ( Figure 1). The climate is semiarid to arid, with average annual precipitation of 70 -100 mm (El-Isa et al., 1995;Salameh et al., 2018).Three   The lacustrine deposits comprise alternating, light-dark, laminated to thin layers of carbonates (aragonite, calcite), quartz and clay minerals (kaolinite), with cm-scale, idiomorphic halite crystals. These are interbedded with 10 localised thin to thick beds of halite and/or gypsum (Figure 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 (Figure 2c and e), which are interbedded with minor silts and clays (El-Isa et al., 1995;Polom et al., 2018;Sawarieh et al., 2000;Taqieddin et al., 2000). At the former stable shoreline, alluvium and lacustrine deposits are in direct lateral contact. 15 Comparable deposits have been widely found along the western shore of the Dead Sea and in numerous boreholes there (Yechieli et al., 1993(Yechieli et al., , 2002. Borehole 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, Figure 1b, Figure 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 20 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 (Figure 3) whose characteristics are summed up in the following 25 based on works of Akawwi et al., 2009;Goode, 2013;IOH andAPC, 1995 andKhalil, 1992. The first is a lower https://doi.org /10.5194/hess-2021-37 Preprint.  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. Figure 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 5 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. p5, Sect. 2.3, IOH and APC, 1995). From these regional aquifers, groundwater water 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 ( Figure 3). The discharge of this spring is around 0.28 m 3 s -1 , and is thought to be principally derived 10 (~ 80 %) from base flow within the Wadi Ibn Hamad (cf. p18, Sect. 5, IOH and APC, 1995), which in turn originates from outcrops of Ram-Kurnub aquifer unit along the wadi bed upstream of Ain Maghara. The remaining (~ 20 %) of 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 (Figure 1b, Figure 2b). Other surface stream channels in the lacustrine deposits lack these 15 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 is not recorded. Head elevations in the vicinity of Ain Maghara are reported to be -300 to -350 m msl in 1994 (cf. pg 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 20 depth was apparently recorded for BH2 (El-Isa et al., 1995). A groundwater well just over 5 km to the 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 -1 from 1960 -2010, whilst over the same time period the measured electrical conductivity had increased by 59 µS yr -1 (US Geological Survey et al., 2013).
The first published report of geophysical investigations in the study area is El-Isa et al. (1995,) which includes 25 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 southwest 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, 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 30 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 Figure 3b. The water levels in the boreholes were combined with spring elevation data to determine a hydraulic gradient of roughly 30 35 m km -1 from east to west, following the surface topography.
A second report, by Sawarieh et al. (2000), includes additional refraction seismic surveys, as well as water chemistry, temperature, pH and electrical conductivity (EC) of springs, wells and sinkholes undertaken in February https://doi.org /10.5194/hess-2021-37 Preprint.  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 the locations of groundwater springs feeding channels CM1-3 and CM6 ( Figure 1b). Groundwater samples were analysed for anion and cation concentrations and total dissolved solids 5 (TDS) in the laboratory at the Ministry of Energy and Mineral Resources of Jordan (MEMR). The results are summarised in Table 1. In general, the measured TDS increased from east to west: groundwater sampled at Ain Maghara could be defined as 'fresh' (EC: < 0.75 mS cm -1 ; TDS: < 500 ppm), whereas water rising from the springs closer to the shore would be categorised as 'brackish' (EC: 1.5 -50 mS cm -1 ; TDS: 500 -35000 ppm), aside from their sample 1 which had abnormally high EC and TDS, being more of a 'brine'. The Na/Cl molar ratios of water 10 samples from these springs were relatively low. The EC readings from the field broadly correlate with the proportions of TDS and Na + and Clions measured in the water samples.  Alrshdan, 2012). GPR data acquired with a SIR-10B system near the mudfactory was 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 20 formation were 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 (> 10000 Ω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 25 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 profiles. Two ERT profiles acquired further north in the study area by Frumkin et al. (2011) also imaged high resistivity anomalies which they interpreted to represent subsurface cavities. https://doi.org/10.5194/hess-2021-37 Preprint. Discussion started: 1 February 2021 c Author(s) 2021. CC BY 4.0 License.
All ERT and GPR surveys on the alluvial fan deposits only penetrated to maximum 25 m deep, 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 -2014. These data were interpreted as incompatible with a proposed 2 -10 m massive salt layer at around 35-40 m depth , but did reveal zones of high S-wave scattering in the 5 subsurface near the main depression as well as potential buried and refilled channel systems within the alluvial fan architecture and a deeper lying "silt & clay" layer at around 60-80 m depth.

Material and Methods
Our multidisciplinary investigation of surface 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 15 ERT, SP, GPR and SH-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 the Appendix A.

Remote sensing and Field surveys
The remote sensing dataset comprises high-resolution optical satellite imagery and aerial photographs that span the 20 years 1967 to 2018. The dataset is identical to that presented in Table 1  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)

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. Although rarely used as a stand-alone method, it is the only geophysical method able to resolve groundwater flow in many circumstances and on a large 5 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 physio-chemical mechanisms are of 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. Corwin, 1990;10 Overbeek, 1952;Revil et al., 1999aRevil et al., , 1999b. See Appendix A for further background information. Our SP survey was conducted in 2018 (for location, see Figure 1b). We used the fixed base configuration ( Figure   4a), with a high impedance (> 50 MΩ) Voltcraft voltmeter and several non-polarizable Ag/AgCl electrodes for the survey (Figure 5e). The semi-permeable bottom filters of the electrodes are covered by a thick, wet bentonite mass 15 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. 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 fresh water 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. 20 More details on the SP array can be found in the results (Sec. 4.2.1).

Electric Resistivity Tomography
Electric resistivity tomography is an active geophysical method within the family of geoelectric methods.
Geoelectrics have been 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.
Nowadays, ERT works by injecting direct current (DC) or low frequency alternating current (AC, short time DC 5 application) into the ground and measuring the differences of the electric potential field at multiple electrodes 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;Jardani et al., 2012;Loke et al., 2018;Loke and Dahlin, 2002). Small electrode spacing leads to high resolution and shallow penetration depths, whereas large electrode spacing lowers resolution and increases 10 penetration depths. More details on these geometric considerations are given in the Appendix A.
Our ERT survey was performed in October 2018 by using a Lippmann 4point light 10W direct current instrument (Grinat et al., 2010) and stainless-steel electrodes connected to the Lippmann multichannel cable system ( Figure   5b-d). Data were collected by using a multichannel control system implemented in the Geotest v. 2.46 software 15 (Rauen, 2016)  coupling. Remaining anomalous data points (showing e.g. physically implausible negative resistivities or nonreliable local high resistivities) have been carefully removed from further analysis.

Shear wave reflection seismics
This method uses the principles of reflection seismology (Sheriff and Geldart, 1995), which is widely applied in hydrocarbon exploration. A seismic source at or close to the Earth's surface generates an elastic wave signal that 30 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 ( Figure 4c). 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 35 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 pore space content (fluids, gases) and they propagate significantly slower compared to compressional body waves (P-waves). The resulting smaller wavelengths of S-waves can https://doi.org/10.5194/hess-2021-37 Preprint. Discussion started: 1 February 2021 c Author(s) 2021. CC BY 4.0 License. improve the subsurface resolution by up to ten 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.
The shear wave seismic surveys reported on here used horizontally polarized S-waves (SH-waves) and were undertaken in October 2014 and December 2016. Details of the data acquisition and processing applied regarding 5 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 perpendicular 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. 10

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. The method typically achieves vertical resolutions between 1-10s of centimetre depending on antenna frequency (Jol, 20 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 general analysis various attributes of karst systems, often in combination with other geophysical methods (e.g. Carbonel et al., 2014;Čeru and Gosar, 2019;Gómez-Ortiz and Martín-Crespo, 2012;Gutiérrez et al., 2011;Kaufmann et al., 2018;Kruse et al., 2006;Mount et al., 2014;Sevil et al., 2017). The method has proven especially useful to determine the limits between unsaturated 25 and saturated zone, image the water table in karst systems, and define the fresh-saltwater boundaries in coastal aquifers (e.g. Gustafsson, 2005;Igel et al., 2013;Kruse et al., 2000;Mount et al., 2015;Paz et al., 2017). https://doi.org/10.5194/hess-2021-37 Preprint. Discussion started: 1 February 2021 c Author(s) 2021. CC BY 4.0 License.
An electromagnetic impulse of a distinct frequency f is sent via a transmitter antenna into the ground, the ground response waves are recorded by a receiver antenna mounted usually close to the sender (so called vertical incidence case, Figure 4d). The survey is taken 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) versus horizontal distance. 5 All GPR surveys were conducted in October 2014 on 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 ( Figure 5f). Several transects were collected using both 100 MHz and 50 MHz unshielded antennas, although only results from the later are reported here for a total of two transects (GPR 1 and GPR 2 in Figure 1b) 10 that followed ERT surveys. Time-to-depth conversion for radargrams was based on estimates of EM wave velocities from three different approaches: 1) from diffractions in common offsets, yielding velocities ranging between 0.09 m ns -1 for deeper materials to 0.15 m ns -1 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-0.2 m ns -1 ; and 3) from samples collected in the field at approximately 2 m deep in a sinkhole wall and measured in the 15 lab under different moisture contents, reaching average values of 0.15 m ns -1 for the driest (5% moisture content) conditions (Appendix A). The mean dielectric permittivity value hence lies around ~ 4, typical for dry sand.

Inversion of SP and ERT data 25
In this study, all SP anomalies were assumed to be generated by a 2D polarized sheet that is inclined, thin and 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 (xo); its depth (h), the inclination angle (αi) and the halfwidth 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 30 algorithm was first described by Eberhart and Kennedy (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. 35 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 https://doi.org/10.5194/hess-2021-37 Preprint. Discussion started: 1 February 2021 c Author(s) 2021. CC BY 4.0 License. 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 5 cells to minimize the difference between calculated and measured apparent resistivity (Loke et al., 2018;Sjödahl et al., 2008).

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 10 Sea level and of the hypothesised 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 2 km long transect perpendicular to the shoreline of the Dead Sea (i.e. along a roughly Norwest-Southeast orientated profile following the NW section of profile AA' in Figure 3, with the centre approximately at the head of canyon CM5). The model salinity distribution provided electrical resistivity values 15 that were compared to those estimated from the ERT inversions. The 2D hydrogeological model setup and approach are presented in Figure 6.
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 25 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;Polom et al., 2018;Taqieddin et al., 2000).
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 hydraulic properties in the model are mostly estimates 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 ℎ = 1 * 10 −6 m s -1 for alluvium and ℎ = 1 * 10 −10 m s -1 for lacustrine sediments are taken as mean values from literature as adapted for the DS in studies from Sachse (2015), Strey (2014) and Shalev et al. (2006). Vertical 10 hydraulic conductivities are assumed to be 10 times less the horizontal ones, to represent anisotropy imparted from sedimentological layering. Hydraulic conductivities for the Dead Sea and for the hypothesised conduit system have been selected after various parameter tests. We added 1 m wide horizontal conduits of high effective porosity and hydraulic conductivity (Table 2), connected by 1 m wide vertical conduits, to the system in each second year of simulation. Due to the groundwater level decline and shoreline recession, the horizontal extent of later-formed and 15 deeper lying conduits is greater that of earlier-formed and shallower conduits. Conduit effective porosity and hydraulic conductivities were first varied in a trial & error approach to achieve model convergence and realistic steady state hydraulic results, as further explained below in Sec. 4.3.
We simulate density driven groundwater flow based on salt concentration in the porewater. The modelling approach 20 is to first simulate the saline water distribution in an estimated steady state over hundred years before the Dead Sea started retreating, and to use this as initial condition for simulating the yearly lake recession. The starting salt concentration of 340 g/ -1 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 https://doi.org/10.5194/hess-2021-37 Preprint. Discussion started: 1 February 2021 c Author(s) 2021. CC BY 4.0 License. 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 water inflow of initial salt concentration 0.67 g/ -1 through the alluvial sediments, adjusted to values of the Ain Maghara spring (Table 1). A hydraulic head gradient of ~ 30 m km -1 (cf. Sec. 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 5 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 1/100, the same holds for changes in the concentration. Formation factors, effective porosities and the salt diffusion coefficient are known for the Dead Sea sediments analyzed samples of the western side of the lake (Ezersky and Frumkin, 2017;Yechieli and Ronen, 1996). Transport steps and dispersivities have been derived from different saltwater intrusion studies and the 10 classical Henry problem (Abarca et al., 2005;Croucher and O'Sullivan, 1995;Geo-Slope, 2004;Langevin and Guo, 2006;Meyer et al., 2019;Zidane et al., 2012). For an overview of the initial spatial distribution of the effective porosity, concentration and head parameters please refer to the Appendix B.  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 versus salinity relationships developed from borehole studies (Ezersky and Frumkin, 2017;Yechieli et al., 2001) and refer to the Total Dissolved Solids (TDS) content rather than the chloride concentration 5 in our simulations.
We can derive the bulk soil resistivities via = ⁄ by using Archie's law (Archie, 1942) even for the claycontaining DS sediments as they do not present cohesion and cation exchange effect after Ezersky and Frumkin 10 (2017).

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 (Figure 7), which is coincident with the area in 15 which the new geophysical surveys were conducted (see Figure 1b for 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/mud-flat boundary or within the mud-flat deposits, and they originated at points that lie over 100 m seaward of -and several meters below -the 1967 lake 20 level. The heads of straight channels do not correspond with spring points, but they have initially developed some distance out on the mud-flat, 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 lakebed as the shoreline retreats over time. While the straight 25 channels also show upstream growth (e.g. CS1-3 in Figure 7), most meandering channels show little or no upstream growth (e.g. CM1-4 and CM6 in Figure 7). Established sections of both channel types also widen progressively  (Figure 7). This depression formed on a scale much larger than the 5 individual sinkholes, and much of its perimeter is delimited by numerous ground cracks and faults . 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 10 groundwater outflow, is CM5. This formed in the summer of 2012 with its head initially located in the middle of   These subsurface channels 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 mudflats (Figure 9a-b). 15

ERT & 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. Figure 8; Al-Halbouni et al., 2017;Watson et al., 2019). These rapid geomorphological changes, occurring at the interface between alluvial cover and lacustrine 5 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 here (Figure 1). The results of both surveys are shown in Figure 10.
Using the PSO 2dD inversion method (cf. Sec. 3.3.1), we extracted and inverted data for five profiles within 10 traversing the SP array (SP1-5, c.f. Figure 10a). The results and fit are shown in Figure 10d. 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 on the location of each anomaly (in x-y space, cf. Figure 10a).
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 15 flow directions, is illustrated in Figure 10e. The length of each strip corresponds to the half-width of the matching anomaly in Figure 10d. The direction of water flow along a strip to produce the inverted anomalies is generally from minus (red) to plus (black). It should be emphasised that the modelled strips are by no means the only possible 'strip' 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 strip. For clarity, the results of profile SP3 have 20 been moved to the Appendix A, which contains also 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 = 50 mV, where = − ( being the potential at the reference electrode).

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 which crosses the boundary between alluvium and lacustrine sediments next to the channel CM1, which was dry in 2018 (although it was weakly active during field campaigns in the years 2014-2016). The results of this survey are shown in Figure 11, which reveals, similar to ERT lines 3 and 4, an 5 extremely low resistivity layer ( = 1 − 5 Ωm) which again coincides with surficial lacustrine mud deposits. The ground above this layer in the south of the profile is more resistive but the spatial distribution of these resistivities 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 (Figure 11d) primarily comprises a light brown marl composed of carbonate and clay minerals, with some laminated horizons 20 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 25 in the intervening pore spaces (cf. Bookman et al., 2004).

10
This degree of cementation is remarkable for the alluvium and is likely responsible for the elevated resistivity of the shallow subsurface at this location.

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;Polom et al., 2018;Sawarieh et al., 2000), and to investigate possible changes to the subsurface hydrology 15 using electrical methods, we performed two new ERT surveys on the alluvial fan, at the southern edge of the uvala (for locations see Figure 1). ERT lines 1 and 2 are coincident with seismic profile lines 1 and 2 as reported by Polom et al., 2018, respectively (Figure 12a-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 southwest

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 Figure 13a, 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 5 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. Figure 12c). 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 (Figure 1b). An extensive area of dense vegetation 10 directly at the former Dead Sea shore-line is nonetheless indicative of a sustained source of sub-surface fresh water 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 ( Figure 13a). The seismic section in Figure 13b

Water table inferred from GPR 10
GPR common offset profiles were processed according to Sec. 3.2.4 prioritizing the enhancement of deep reflectors in order infer the location of the water table during data collection in Oct. 2014 (Figure 14a and b). Both profiles show the presence of shallow (0.3 m ns -1 ) 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, Figure 12c and d). GPR profile 1 shows two major air 15 reflections (i.e. hyperbolas, black arrows in Figure 14a Figure 14b as the water table. It appears more horizontal than in GPR 1, which may indicate the absence of strong lateral variations. When compared to the first 280 m of ERT 2, this depth corresponds to the transition zone between moderate and low resistivity (10 -100 Ωm). 15 and dry (in situ) alluvium ( ~ 5 %), respectively (53Appendix A).

Hydrogeological modelling
To study the effects of the falling Dead Sea level upon the subsurface hydrogeology in the study area, we applied  The main criteria for evaluating the quality of the hydrogeological model solutions were: (1) model convergence; (2) realistic hydraulic head development without singularities or jumps, which is expressed as a smooth water table After the initiation period simulation of 100 years under steady state conditions, the salinity of the system ( Figure   15b) is initially in a quasi-equilibrium with a clear distinction of fresh and saline water according to the local geology. The fresh-saline water system reacts relatively slowly to base level fall because the diffusion of the saline water is a slow process. As such, the reaction of the hydrogeological system is visible only after 15 years (or 15 m 35 vertical lowering of the DS). The acceleration is caused by the development of a network of conduits. Formerly saturated areas become dry, and a rather salty water composition dominates the flow in the conduit system with concentrations between 125 and 250 g l -1 . The fresh-saline interface system builds locally a complex shape, and https://doi.org/10.5194/hess-2021-37 Preprint. Discussion started: 1 February 2021 c Author(s) 2021. CC BY 4.0 License. different levels of saline water can be encountered in a hypothetical vertical slice, but also lateral variations appear.
Regionally it evolves into a landward inclined shape of a Ghyben-Herzberg style interface.
The soil resistivity shows a duality of responses as the model develops, increasing in the near-surface lacustrine sediments (due to renewed inflow of fresher groundwater in the conduits), but decreasing gradually (Figure 15c) 5 at the diffusion fronts. The latter is especially true for the alluvium-mud boundary, but also in the implemented channel water system ranging far inwards the alluvial sediments. Due to saline water intrusion inland, the resistivity is gradually lowering in the deeper areas of the model, including areas of alluvium that were initially saturated by fresh water (US Geological Survey et al., 2013).

15
At the contact between the alluvial and lacustrine sediments (at a distance of 1500 and 1900 meters from the shore; For comparison with resistivity structures observed in our ERT surveys, Figure 16b shows the slice of the simulated bulk soil resistivity results for saturated areas. We can see a typical evolution at the contact zone between very 25 conductive lower layers and very resistive upper layers. Once the saline water intrusion has developed further, the resistivities are modulated gradually depending on salt concentration and material type. The conduits appear usually as more conductive lines in both material types.

Discussion
Here we first discuss the geophysical results and their interpretation, and we then discuss 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 at Ghor Al-Haditha sinkhole area. 5

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. 10 Small variations in the portions of metal (e.g. Fe) or salt (e.g. Na, K), whether contained within in the solid or liquid phases, can result in large variations in their electrical conductivity. Additionally one has to bear in mind the 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 as a thicker, less conductive layer. The conductance of both features would image similarly, and interpretation 15 needs to be done carefully, on the basis of additional data, boreholes and conductivity measurements.
Consistent with former studies in the area (Alrshdan, 2012;Frumkin et al., 2011), our geophysical surveys on the alluvial fan (ERT1 and ERT2, Figure 12) indicate a general decrease in resistivity with depth, broadly stratified into distinct regions. On the basis of logs of sinkhole walls nearby (Sawarieh et al., 2000;Taqieddin et al., 2000), 20 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 (Figure 12c), which may represent initial sediment geometry, typical different saturation grades and/or salinization (e.g. Farzamian et al., 2019aFarzamian et al., , 2019bGonçalves et al., 2017). This 25 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 30 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 35 and 0.4 Ωm for lacustrine sediments, and between 0.45 and 5.8 Ωm for alluvial fan sediments (Ezersky and Frumkin, 2017;Yechieli et al., 2001). 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 https://doi.org/10.5194/hess-2021-37 Preprint. Discussion started: 1 February 2021 c Author(s) 2021. CC BY 4.0 License. parts (with resistivities lower than 5 Ω ) 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 (Figure 12c) 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 5 transition between resistive and conductive areas in the profiles and in the models offer indications of the freshsaline 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. The GPR profiles 1 and 2 10 hereby help to constrain the limit between saturated and unsaturated ground of 30-38 m below surface in October 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 25 road on profile ERT1 (Figure 12c; Sec. 4.2.4). Others are found in association with the SP anomalies upslope of the stream channel head, near the interpreted base on the alluvium and just above the underlying lacustrine deposits (Figure 10b-c). The suggestion that these low-resistivity ERT anomalies represent flow in such a conduit network is consistent with electrical conductivity measurements of the springs at the head of CM5 during the 2018 field campaign. These measurements, between 22 -26 mS cm -1 , would equate to resistivities between 1 -0.5 Ωm. 30 Similar resistivities were measured both at the low-resistivity anomaly at the edge of the uvala, and at around 420 mbsl 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.

35
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.
can develop locally at the alluvium/mud contact, as indicated by the SP results ( Figure 10) and by the simulated flow pattern (Figure 16a). 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).

Feasibility of the hydrogeological model
The nature of aquifers formed at a boundary between fresh and saline water is critically dependent on the physical 5 and chemical contrast between saline and fresh water. Fluid flow driven by density contrasts between saline and fresh water is generally described by two mathematical approaches: (1) the 'sharp-interface approximation' (Ghyben, 1888;Herzberg, 1901;Pool and Carrera, 2011), which assumes two fluids of constant density, no mixing and a pronounced saline-fresh water interface; and (2) the density-driven approach, which accounts for mixing of the two water types within a 'transition zone' (Croucher and O'Sullivan, 1995;Henry, 1964;Simpson and Clement, 10 2004). Hydrogeological modelling of groundwater flow using both approaches has been previously undertaken for sites around the Dead Sea area (Alfaro et al., 2017;Odeh et al., 2015;Sachse, 2015;Salameh and El-Naser, 2000;Shalev et al., 2006;Strey, 2014;Yechieli, 2000), 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 by Strey (2014), for the Darga Quaternary Dead Sea group 15 sediments, with a considerable analysis of factors controlling model performance, and Shalev et al. (2006), with an application to salt dissolution and sinkhole formation by faults as preferred groundwater conduits.
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 20 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. 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 25 sedimentation, evaporation and erosion periods (e.g. Bartov, 2002;Levy et al., 2020;Neugebauer et al., 2015). 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 30 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 35 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, although in homogeneous systems groundwater flow can be assumed perpendicular to the shoreline (effectively 2D), in https://doi.org/10.5194/hess-2021-37 Preprint. Discussion started: 1 February 2021 c Author(s) 2021. CC BY 4.0 License. complex 3D aquifers, groundwater flow directions are deviated by geological heterogeneities and hence angular (Meyer et al., 2018b(Meyer et al., , 2018a. 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 5 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 the 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 10 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 (Bakalowicz, 2015;Bakalowicz et al., 2008;Fleury et al., 2007). At Ghor Al-Haditha, our surface manifestations indicate that such conduits can be generated by physical erosion of the weak alluvial and lacustrine materials, by chemical erosion of the evaporite component of the lacustrine deposits 15 or by both processes. These processes are not explicitly simulated in the model, but incorporation of these and the related mechanical feedbacks could be 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 20 (Ezersky andFrumkin, 2017, Yechieli, 2006). The approach is based on the definition of empirical parameters, the formation factors, a classical issue in hydrogeology. 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, 25 erosion would enlarge these conduits, such that they may become more electrically conductive or resistive,

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 (Arav et al., 2020;Avni et al., 2016;Shviro et al., 2017). 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). 5 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 thus promotes the local development of sinkholes, which in turn feed surface water down into evaporite karst conduits, which undergo chemical erosion and destabilise and form more sinkholes in a self-accelerating manner. This 10 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 15 of individual flood events as at Ze'elim, 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 southwest to the north-east of the study area over the past 30 years with a focussing of groundwater discharge at CM5 from 2012 onward (Figure 1b, Figure 7). This broadly coincident with a shift in the area of active subsidence from the 20 south-west to the north east also . 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 southwest. 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 that area has been minimal. Only since 2006 has surface run-off from the Wadi Mutayl 25 shown signs of being partly directed into the uvala (Figure 7), 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 Wadi Mutayl into the uvala may have contributed to its accelerated development during (cf. Avni et al., 2016.
Nonetheless, the observed shifts in discharge location and subsidence around the main Ghor Al-Haditha alluvial fan must overall reflect a re-routing of groundwater flow within the alluvial deposits and any potentially underlying 30 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. Brook and Murphy, 2017;Despain and Stock, 2005;Farrant and Simms, 2011;Palmer, 1975;Plan et al., 2009;Simms and Hunt, 2007;Šušteršič, 2006). New conduits, often termed 35 '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 https://doi.org /10.5194/hess-2021-37 Preprint.  of the study area and is slowest in the vicinity of CM5, relative to the 1967 shoreline (cf. Figure 3 of Watson et al., 2019).
Our study also complements the recent works 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 5 event of 2015, discharge at the main spring (s4) clearly emanated from a conduit within the lacustrine sediments, about 6 m below the alluvium/lacustrine interface (Figure 9e-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 (Arav et al., 2020;Avni et al., 2016), high flow events were inferred to cause substantial erosion and rearrangement of the conduit network. This erosion seems to be largely chemical (i.e. 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 earlier studies at the Dead Sea (Al-Halbouni et al., 2017;Arkin and Gilat, 2000;Taqieddin et 20 al., 2000). 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 base flow following recharge of the fan by rainfall and/or floods, and they can be sustained by the cementation of the alluviumespecially the older deposits. It is stressed that such a physical mechanism of conduit generation 25 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. 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.

30
The spatial extent and dimensions of any conduit network at Ghor Al-Haditha remains 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 in Watson et al. (2019), we regard dissolution by 35 surface water to be almost absent at Ghor Al-Haditha. In our study area, the medium which provides the initial hydraulic connectivity are the interpore spaces in the matrix. 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 (Klimchouk et al., 2000;Palmer, 1991Palmer, , 2007. No direct https://doi.org /10.5194/hess-2021-37 Preprint.  observations of spongework conduit zones have been made in evaporite karst areas and they are rare in the case of limestone karst. The locations where spongework has been observed are often areas of mixing of different water chemistries in rocks with high primary porosity, such as mixing of fresh water and seawater along the coastal carbonate platforms of the Bahamas, where they are termed 'flank-margin caves' (Mylroie and Carew, 1990); the coasts of the Yucatan Peninsula in Mexico (Back et al., 1984;Smart et al., 2006); and the Prichernomorsky Basin 5 on the north coast of the Black Sea (Klimchouk et al., 2012). In the case of the Prichernomorsky Basin, water derived from deep sources rich in sulphuric acid mixes with shallower matrix flow, causing aggressive groundwater erosion under hypogenic conditions. We suspect that deep-seated groundwater sources rich in sulphur 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 s2 at CM5 (Figure 8) but particularly those to the north of the studied 10 area such as that feeding CM7, have a highly sulphurous odour. It has been noted that other groundwater springs in the vicinity smelling of sulphur seem to derive from the Ram and Kurnub sandstone units (IOH and APC, 1995;Khalil, 1992), which would likewise imply a deep-seated source to these waters at Ghor Al-Haditha.

Conceptual model of hydrological and geomorphological processes
The culmination of our study is a refined conceptual model of the surface and subsurface hydrological processes 15 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 ( Figure 17). The fall in base level in the study area has led to the development of a strong hydraulic gradient and an ingress 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 within the smaller Wadi Mutayl. Leakages from the bedrock 20 aquifers may also play a role, especially in the northern part of the Ghor Al-Haditha sinkhole site (north of the area shown in Figure 1b), where the prevalence of sulphurous 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 deep these 25 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 (Figure 13).
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 30 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 (Figure 12c) 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-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 35 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 https://doi.org/10.5194/hess-2021-37 Preprint. Discussion started: 1 February 2021 c Author(s) 2021. CC BY 4.0 License. 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 northeasterly shift over time (Figure 1b, Figure 7), 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 5 levels. The shallow subsurface is dominated by 'advection' transport processes within the matrix and conduit network, while the late initiation of saline water intrusion only after 17 years (or 850 m) of DS regression in our model (bottom row of Figure 15, Sec. 4.3) suggests 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 Ω may extend remarkably more towards 10 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 more comprehensive information on the subsurface lithology 15 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 to 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 ground truth of the geophysical surveys undertaken in the study area and to improve our understanding of the nature of erosion in the subsurface.   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). Line X-Y is the location of the cross-section in part (b), which shows the 5 hydraulic structure of the subsurface as inferred from our geophysical studies and the hydrogeological modelling.

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 -30 m deep, bringing fresher groundwater to 10 the constantly-retreating Dead Sea shoreline. 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.

Data availability:
All geophysical data are available on request from the authors. Free satellite images (Corona) and photogrammetric 25 survey raw images, DSMs and orthophotos are available upon consultation with the authors. Geological Map 1:50,000 Ar Rabba: 675 available at discretion of MEMR.

Video supplement:
Three videos of the karstic spring s4 formation at canyon CM5 and associated upstream canyon growth by 30 retrograde erosion are provided in the online supplement.

Author contributions:
DAH led the conception and writing of the manuscript. DAH, RAW and EPH designed and planned the field experiment and largely assembled the manuscript. DAH and RAW created and assembled the figures. DAH 35 analysed the ERT, SP and GPR data; FDS analysed the SP data; RAW analysed RS data; UP analysed the seismic data; XC analysed the GPR data. RM supervised and verified the numerical simulations. RAW and HAR participated in field data collection. DAH conceived the original idea. CMK and TD supervised the project. All authors contributed to improvement and editing of the manuscript.  We here present details on the different geophysical methods and additional supportive results.

Self-Potential:
The following background on Self-Potential is based mainly on comprehensive works from (Jardani et al., 2007;Jouniaux and Bordes, 2012;Jouniaux and Ishido, 2012;Richards et al., 2010;Vichabian and Morgan, 2002). For 15 more details, the reader is referred to these authors and references within.
Clay minerals carry negative superficial electrical charges due to crystallographic defects. Attracted positive charges (cations) stay therefore 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 20 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 negative zeta-potential for the movement of positive charges and forms the physical 25 background for the calculation of the self-or streaming potential. We derive the governing main equations for streaming potential calculation in the following. The coupled flow equations for saturated case are after (Sill, 1983): Jf is hereby the primary flow of matter (e.g. Darcy flow of water) and the effects of the coupled secondary flow Je are usually considered as negligible. However, here Je 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 30 and in the grain. Lek is the electrokinetic coupling coefficient following Onsager´s reciprocal relation in steady state conditions. ηf is the dynamic viscosity of the fluid, σ0 its' specific conductivity, k0 the hydraulic conductivity and ∇P the pressure gradient. 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: C = ∂V/∂P = εf ζ/ηf σeff dominates. Typical SP gradients of groundwater flux along an axis x are around dV/dx = 0.2 mV m -1 . For karst aquifers this equation holds for any capillary and the SP signal comes from capillaries and 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 e.g. sinkhole downward seepage.

10
The detailed results of the inversion of five SP profiles extracted from the SP array are summarized in the following table, the result for SP3 is presented in Figure A1. 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 structure. The dipole moment for inversion has been constrained to positive values only and relates the current density J to resistivity for a thin 15 polarized strip like structure to = 2 ⁄ . For more details on SP modelling and inversion see (Biswas, 2017;Monteiro Dos Santos, 2010;Murthy and Haricharan, 1985).  Estimation of the penetration depth however is difficult subject to the concept of conductance vs. thickness of a 10 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 manuscript (Sec. 3.3.1).
Ground penetrating radar and sample laboratory tests:

15
Water content may lead to strong damping of the electromagnetic signals by direct current losses, preventing sufficient depth penetration of GPR signals. Therefore, laboratory measurements have been performed to determine the properties of alluvial material samples from the field site. Figure A2 contains the results for electromagnetic laboratory tests on the samples with a frequency of = 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 20 area, the conductivity = 1/ lies between 2.7 mS m -1 and 0.5 S m -1 . The relative dielectric permittivity ε lies between 4 and 24, depending also on the moisture content. This results in an attenuation of 10-100 dB m -1 for 100 MHz and the penetration depth δ, or skin depth, for a 100 MHz antenna in a homogeneous half space can be approximated to δ = √ 2 (Jol, 2009)

5
The frequency 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 one has to accept the low horizontal resolution of Δ = √ 2 with r as the radial distance from the sender. For 50 MHz and 0.15 m ns -1 velocity the minimum horizontal extension of an object that can be resolved in 10 m depth is 3.87 m. For 100 MHz it is 2.7 m. 10 Appendix B

Hydrogeological model parameters, alternative models and resistivity calculation 15
The hydrogeological modelling procedure as described in Sec. 3.3 requires the definition of further important parameters such as effective porosity and concentration distribution, initial conditions and fixed head and concentration locations ( Figure B1).
Alternative models for a convergent steady state solution have been tried and are presented in Figure B2. There are 20 no channels installed in these models, rather only the hydraulic conductivities of lacustrine sediments are varied. Note that in model 1 the water table does not follow a realistic curve. In model 2 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 one order of magnitude higher than AS material, an unrealistic scenario for the whole lacustrine material. Hence, these models were depreciated in favour of the model presented in Sec. 4.3 with local channels of high permeability.